Comparison of density functional and MP2 geometry optimizations of Na(H2O)n (n = 1–3) clusters

THEO CHEM ELSEVIER

Journal of Molecular Structure (Theochem) 425 (1998) 87-94

Comparison of density fknctional and MP2 geometry optimizations Na(HzO), (n = l-3) clusters

of

Daniel Bacelo, Yasuyuki Ishikawa* Department

ofChemistry, University of Puerto Rico, PO Box 23346, San Juan. PR 00931-3346. US.4 Received 2 November

1996; revised 4 March 1997

Abstract Geometry optimizations have been carried out for clusters of sodium with as many as three water molecules. Potential energy surfaces were generated both conventionally, via Hartree-Fock calculations with electron correlation correction by secondorder Merller-Plesset (MP2) perturbation theory, and by three methods based on density functional theory (LIFT). The DFT methods have proven to be able to reproduce the most stable geometric configurations of the hydrogen-bonded clusters predicted by MP2. Some discrepancies occurred with the DFT methods, including erroneous energy ordering of local minimum structures and location of false minima on some potential surfaces. 0 1998 Elsevier Science B.V. Keywords:

Density functional theory; Potential energy surface; Ab initio correlated calculation; Sodium-water

1. Introduction In the last decade there has been research interest in the structure, properties and dynamics of solvated species, and in the spectroscopy of solvation in bulk liquids and clusters [ 11, interest spurred by the importance of these systems in chemical and biological processes. In particular, studies of solvated alkali metal atoms [2-51 by photoionization mass spectrometry have focused attention on the need to be able to predict the properties of these systems theoretically. Hydrated sodium atom clusters have served as prototypes. Curtiss et al. [6] studied the fundamental structure and bonding characteristics of the NaH20 complex. Hashimoto et al. [7] examined sodium solvated by as many as six water molecules in ab initio unrestricted Hartree-Fock calculations, and * Corresponding Ol66-1280/98/$19.00 PII

author. 0 1998 Elsevier

SO166-1280(97)00132-2

Science

clusters

Bamett and Landman [8] have studied a similar range of these complexes using density functional theory (DFT). These last two studies both focused on finding cluster ionization potentials in order to compare theoretical predictions with experimental mass spectrometric results. Such determinations require the cluster energy at the minimum energy configuration, and the studies thus sought those configurations. They also reported finding several metastable configurations, local minima on the cluster potential surfaces. Theoretical determinations of the structures of stable cluster isomers are valuable because such structures are not experimentally accessible, yet they are necessary for interpreting cluster dynamics. Recently, Ishikawa et al. [9] attempted to determine all the minima, global and local, for sodium clustered with as many as three waters, using ab initio Monte Carlo simulated annealing. The structures which were located were then energy optimized in

B.V. All rights reserved

88

D. Bacelo, Y. IshikawdJournaI

of’blolecular

conventional basis set Hartree-Fock (HF) calculations with electron correlation correction by secondorder Moller-Plesset (MP2) [lo] perturbation theory. Molecular geometry optimization at the MP2 level of theory is widely employed. Geometries and harmonic vibrational frequencies derived by gradient analysis of MP2 potential surfaces are reliable [ 11,121. In addition, its deficiencies are understood, and techniques for incremental improvement of MP2 results are available. Methods based on density functional theory have emerged as useful tools for calculating the molecular structures of large systems, especially in the study of metal clusters [ 131. DFT methods are attractive because they offer savings in computational time over conventional HF plus MP2 calculations [14], particularly when the local spin density approximation (LDA) to exchange and/or correlation is used. Progress in recent years in the development of nonlocal corrections to the local density approximation [15- 171 has produced DFT methods able to predict molecular geometries reliably. Hydrogen bonding plays an important role in the stability of hydrated sodium clusters [7-91. The number of distinct isomers and their relative stability depend on hydrogen bonds among water molecules. Because hydrogen bonds are weak relative to chemical bonds, cluster isomers are close in energy, and many are likely to be thermally accessible at normal temperatures. For the same reason the quantitative accuracy with which sodium-water clusters are described is sensitive to the level of theory and to the basis sets used. The various DFT electronic structure methods differ in the forms of the functionals used to describe electron exchange and correlation effects, major electronic effects. The ability of some DFT approximations to describe weakly interacting systems has therefore been questioned [ 181. However, Novoa and Sosa [ 191 have shown recently that the B-LYP and B3-LYP methods [20,21] based on the exchange potentials of Becke [ 16,201 and the Lee, Yang, Parr [ 171 correlation potential do predict the structures of small hydrogen bonded dimers well. Their observation offers promise that DFT methods able to accurately treat systems with extensive hydrogen bonding can be devised. In the present work we extend a previous MP2 study [9] to include an examination of the performance of

Structure (Theochem) 425 (I 99%) 87-94

DFT functionals in predicting the structures and energetics of hydrated sodium. The object of the research has been to determine how well methods based on both local exchange-correlation and nonlocally corrected DFT functionals reproduce the cluster geometries and energies predicted by MP2 optimization, for a given basis set size.

2. Computational

methods

The calculations were done with the GAUSSIAN 94 [22] system. 6-31G* [23] basis sets were employed, with the expectation that they could produce quantitatively accurate bond lengths and angles. Honegger and Leutwyler [24], in a study of water clusters, showed that 6-31G* basis sets do predict the geometries of those clusters accurately at the HF level. MP2 “frozen core” geometry optimization calculations for Na(HzO), (n = l-3) have been previously reported [9]. At the MP2 minimum energy configurations of those calculations, single-point frozen core MP4 calculations with larger 6-31 lG** basis sets were performed to confirm the energy ordering of the stable species. Additional MP2 exploration of the potential surface of Na(H20)3 was done in the course of this study. DFT calculations were performed with three approximations: exchange-correlation the LDA Slater-Volsko, Wilk, Nusair (S-VWN), and the gradient-corrected B-LYP and B3-LYP approximations. The S-VWN potential consists of the Slater local exchange expression [25] with the parametrization by Vosko, Wilk and Nusair [26] for the correlation functional. The B-LYP approximation incorporates gradient corrections into the exchange and correlation energy approximations. The gradient, or nonlocal, correction to the exchange energy suggested by Becke [ 161 is used, together with the correlation functional proposed by Lee, Yang and Parr [17]. Calculations with one of the recently suggested hybrid HF-DFT methods, the B3-LYP [20], were also performed. In the B3-LYP potential the exchange potential is a weighted average of Slater and HartreeFock terms together with Becke’s [ 161 nonlocal exchange. Correlation is provided by the functional of Lee, Yang and Parr [ 171.

89

D. Bacelo, Y. Ishikawa/Journal of Molecular Structure (Theochem) 425 (1998) 87-94

Table 1 Optimized

geometries

of the NaHzO clustera MP2

Parameter

W&O)

0.9741 0.9741 2.3622 104.9373 115.4275 115.4489 128.2498

R(H,-0) R(Na-0) A(H,-OpHb) A(H,-0-Na) A(HbpO-Na) D(NaaO-Hh-Ha) ’ Distances in Angstroms, in Fig. I.

B-LYP

B3-LYP 0.9760 0.9761 2.2866 104.7096 112.2811 111.9575 12 1.9040

angles in degrees. All optimizations

0.9895 0.9896 2.2965 103.7566 109.2871 108.7632 116.2514

were performed

LDA 0.9872 0.9873 2.2129 104.3864 107.8426 107.4387 114.3390

with 6-3 1Ci’ basis sets. Atomic subscripts correspond

to those

3. Results and discussion

3.1. NaHzO

The significant geometric parameters of each cluster structure of one sodium atom with from one to three waters are shown in Tables l-3, respectively, while Figs. l-3 display the structures. In general, the bonds, that is chemical bonds and particularly hydrogen bonds, are predicted to be slightly shorter by the DFT methods than by MP2 with the same basis sets. Thus the DFT methods seem to overestimate the hydrogen bond interactions. MP2 predicts 0-Na-0 bond angles more closed, by typically 3-9, than the angles predicted by any of the DFT approximations. The DFT potential surfaces for Na(H20)2 were found to be more complex than the MP2 in the sense that more DFT local minima were found on two of the surfaces than on the MP2. We discuss the results for each cluster in turn.

Every method used predicted the most stable structure to be one in which the pair of hydrogens is tilted out of the plane defined by Na, 0, and equal projections of the H atoms by angles of 50 to 65”, as shown in Fig. 1. Table 1 summarizes the optimum NaH20 geometries found by each method. The MP2 and DFT results agree well. In particular, the B3-LYP and MP2 results are in excellent agreement. All three DFT methods underestimate the Na-0 distance and overestimate the O-H distances in comparison with the MP2 predictions, but all predict an H-O-H angle in close agreement with the MP2 and the experimental results. Curtiss et al. [6] examined NaH20 in MP2 calculations with basis sets similar to those employed in this study. They carefully examined the nature of the bonding in the complex with electron density and

Table 2 Optimized

geometries

of the Na(H,0)2

Parameter R(Naa0,) R(Na-Ot,) R(H,-0,) R(D,-Hb) R(D,-H,) A(O,-NaaOb) A(Hb-Da-H,) D(D,-H,-Db-Na) d Distances in Angstroms, in Fig. 2(a).

MP2 2.5075 2.3696 2.0886 0.9776 0.9765 71.15 104.1531 7.1977

clustera B3-LYP

B-LYP

2.3793 2.3152 2.0928 0.9807 0.9792 74.1023 104.3996 3.956

angles in degrees. All optimizations

were performed

2.3717 2.3337 2.1327 0.9943 0.9927 75.5206 103.8539 2.2192

LDA 2.3256 2.2389 1.8013 0.9946 0.9913 70.88 105.1476 2.4606

with 6-3 1G’ basis sets. Atomic subscripts correspond

to those

90

Table 3 Optimized

D. Bacelo, Y. Ishikawa/Journal of Molecular Structure (Theochem) 425 (1998) 87-94

geometries

of the Na(H*O),

clusters”

Isomer

Method

R(Na-0,)

R(Na-O,,)

R(O,-H,)

R(O,-Hb)

A@-Na-Ob)

a

MP2 B3-LYP B-LYP LDA MP2 B3-LYP B-LYP LDA MP2 B3-LYP B-LYP LDA

2.4680 2.3579 2.3456 2.2567 2.5217 2.4158 2.4064 2.3555 2.4793 2.3461 2.3420 2.2582

2.3159 2.2636 2.2800 2.1817 2.5215 2.4162 2.4058 2.3508 2.3621 2.2582 2.2778 2.1763

1.9313 1.8898 I .8972 1.6914 2.042 2.073 1 2.144 1.759 1.9714 1.9092

1.7859

96.5212 100.7447 103.3164 103.0872 66.5608 70.0884 72.31 13 65.6376 104.147 101.2638 104.3753 102.0605

b

C

’ Distances in Angstroms, in Fig. 3.

angles in degrees. All optimizations

1.9203 1.7019

were performed

1.7437 I .7467 1.5778 2.0395 2.0684 2.1388 1.7617 1.9551 1.7641 1.7646 1.5922

with 6-3 lG* basis sets. Atomic subscripts correspond

to those

density difference maps and by Mulliken population analysis. They found the out-of-plane angle to be 46.2” and the Na-0 distance 2.343 A, in substantial agreement with the results shown in Table 1. The nonplanarity of the complex offers a sensitive test of theoretical methods. At the MP2 level the complex is nonplanar by only 0.3 kcal mol-‘; Curtiss et al. [6] found 0.24 kcal mol-‘. Much of the energetic difference in the two forms is due to correlation energy. At the HF level the nonplanar structure is favored by only calories. The DFT calculations all favor the nonplanar complex, and all by less than 1 kcal mol-‘. The Na-0 bond in NaH20 is rather weak, 2535% longer than in NaOH, NazO or NaO. Curtiss et al. [6] determined the dissociation energy of the complex to be correspondingly small, the bonding a complex mix of electrostatic, charge transfer, repulsive and dispersive effects. Nevertheless, the DFT methods predict the Na-0 distance well, and they seem to reproduce the subtle effects needed to describe the structure of

the complex. Curtiss et al. [6] found, as expected, significant charge transfer from Hz0 to Na. The Na-Hz0 interaction occurs at the lone pair side of the oxygen, the lone pair pushing charge at the metal atom while itself being drawn toward it. If

Fig. 1. Structure of NaH*O. The Na-0 as connecting rods.

Fig. 2. Structure of Na(H,O)*. The Na-0 and O-H bonds are shown as connecting rods, the hydrogen bond by a broken line.

and O-H

bonds are shown

a H

a

b

D. Bacelo, Y. ishikawa/Journal

Fig. 3. Structural

isomers of Na(H20),.

The Na-0

of Molecular

and O-H

lone pair charge donation is assumed to be a major contributor to the stability of the metal-water bond, the nonplanar form of the complex may be rationalized as the structure which best directs an oxygen lone pair at sodium. The optimum angle of nonplanarity would thereby be predicted to be near 55”, or half the

91

S~~cture (Themhem) 425 (1998) 87-94

bonds are shown as connecting

rods. hydrogen

bonds by broken lines

tetrahedral angle. The calculated values do cluster about this value. At the MP2 level the angle is 52”; angles of 58-67” were obtained with the DFT methods. Note that the water molecule is oriented toward Na in the nonplanar structure as one water is oriented toward the O-H bond of the second in the water dimer.

92

D. Bacelo, Y. Ishikawa/Joutnal

of Molecular Structure (Theochem) 42.5 (1998) 87-94

3.2. Na(HtO)J A structure, displayed in Fig. 2(a), with two Na-0 bonds and a hydrogen bond between the two water molecules, was found to be most stable by all of the theoretical methods employed. DFT methods have been reported at times to have difficulty in representing hydrogen-bonded structures. Barnett and Landman [8] observed that their LDA calculations on Na(H20)* yielded a minimum energy structure which was much too open, i.e. nearly linear, without a hydrogen bond. With gradient exchange and correlation corrections, they obtained a structure near the minimum energy structures obtained in this study. The S-VWN potential is an LDA; nevertheless, it predicts a hydrogen-bonded complex as the most stable configuration. Each DFT method does locate a local structure which is open (see Fig. 2(b)). Such a structure is a saddle point on the MP2 potential surface and on the B3-LYP; harmonic frequency analysis of each surface finds a single imaginary frequency. Frequency analysis of the B-LYP and S-VWN open structures, on the other hand, yields all-positive frequencies, indicating that in these models the no-hydrogenbond structure is a local minimum. The B-LYP method comes nearer to predicting the open, nearly linear structure to be the global minimum; on its surface the open structure lies less than 1 kcal mol-’ above the hydrogen-bonded structure. In addition, the B-LYP surface contains a second open structure, quite similar to the one we have just described, but with the hydrogens of one water molecule flipped approximately trans to their position in the local

Table 4 Energies relative to isomer a (Fig. 3) of Na(HrO)r

minimum structure. This last structure is a saddle point, however, rather than a minimum. The introduction of spurious minima into the potential surface by the DFT methods is most likely an exchange energetic effect, although we have not attempted to analyze the effect in this study. Exchange energies are typically larger than correlation energies, and two of the DFT methods which differed in the number of extra minima produced, B3-LYP with none and B-LYP with one, employ the same correlation functional. The potentials employed are parametrized to reproduce sets of data [ 16,201 which are unlikely to be extensive enough to allow DFT methods to mimic every subtle effect which can arise in the course of theoretical explorations. Of the DFT methods,, the B3-LYP potential did not introduce minima which were not on the MP2 surface. The exchange potential in the B3-LYP potential incorporates HF exchange and was derived to improve DFT description of molecular bond formation and dissociation. In this case the goal appears to have been attained. Becke [20] has pointed out that there is no known systematic procedure for improving radient corrections to the LDA, and that physical intuition forms a large component of what progress is made, because the physics of why they work is not yet fully understood. Table 2 presents a selection of significant geometrical parameters as predicted by each method. Again the B3-LYP results agree slightly better with the MP2 than do the results ofthe other two methods. The hydrogen bond length is slightly overestimated in the BLYP approximation, while it is underestimated with SVWN. There is at least 0.1 A difference in the MP2 and DFT predictions of the metal-oxygen distances.

clustersa

Cluster

MP2 b

B3-LYPb

B-LYPb

LDAb

MP4C

a

0.0 0.3 185 0.8219 2.2443 2.9430 6.9635

0.0 1.953 1 0.5289 1.5514 2.9555 2.2599

0.0 2.4022 0.4813 0.6571 2.6870 0.8618

0.0 0.6817 0.5538 4.7475 3.4450 7.8702

0.0 0.7894 0.93 12

b

c d e f

’ Energies in kcal mol-‘. b Energies calculated with 6-31G’ basis sets. ’ Energies calculated with 6-3 11 G” basis sets

1.7570 2.3864 4.9743

D. Bacelo. Y. lshikawa/Journal

ofMolecularStructure

3.3. Na(H?O)j Six stable isomers were identified for sodium plus three waters [9]. Their structures are displayed in Fig. 3, and Table 4 sets out their relative stabilities as determined by each of the theoretical methods employed. In general, the DFT methods predict larger angles and shorter distances than are found with MP2. Structure 3a was found to be the global minimum by all four theoretical methods, but there are differences in predictions among them concerning the order of stability of the rest of the isomers. Energy differences among the isomers are small, a few kcal mol-’ at most. In order to test the order determined at the MP2 level, single point MP4 calculations were carried out at each minimum with larger 6-3 1 lG*’ basis sets. The MP4 calculations confirmed the MP2 energy ordering. Their effect was to enhance the predicted stability of each isomer. A set of single point calculations with the B3-LYP functional was also carried out with the 6-3 1 lG** basis sets to test the stability of DFT energy ordering with respect to basis set. The calculated order was maintained, and the calculated energy separations changed very little, in every case less than 0.5 kcal mol-‘. During the course of this study, MP2 calculations in addition to those previously reported [9] were undertaken to further explore the potential surface near each minimum structure. An additional local minimum structure was found. The structure is a variant of 3e in which one of the waters which participates in hydrogen bonding is rotated so that the H-bond hydrogen is converted into a dangling one, and the formerly dangling hydrogen assumes the bonding position. The energy difference between these two structures is quite small, about 0.2 kcal mol-‘. The discovery of this subtle axial-equatorial sort of isomerism raises the likelihood that such structures are to be found wherever one hydrogen atom of a water molecule participates in a hydrogen bond while the second does not. Further study is in progress. On the MP2/6-3 1G* potential surface the most stable isomer (3a) is a cyclic structure in which the sodium and the three oxygens form a nearly planar four-membered ring, and the water which contains the oxygen not bonded to Na is held by hydrogen bonds to the other two waters. Structure 3a is only 0.3 kcal molK’ more stable than 3b, 0.8 kcal mol-’ at the

(Theochrm) 425 (1998) 87-94

93

MP4 level. In 3b sodium sits atop a cyclic structure formed by the three waters hydrogen bonded to each other. The most stable structure of a water trimer is just such a cyclic structure as is found in isomer 3b [24]. The doubly hydrogen-bonded isomer 3c is also quite close in energy to the most stable structure, lying only 0.8 kcal mol-’ above it, 0.9 kcal mol-’ at the MP4 level. Isomers 3b and 3c are structurally very similar, 3b with three hydrogen bonds and 3c with two. The position of Na is approximately the same in the two structures, sitting above an O-O-O triangle, equilateral in 3b and isosceles in 3c. 3c has given up some stability with the breaking of one of the hydrogen bonds of 3b, gaining relief of the ring strain of the cyclic trimer. Because the structures are so similar, the two species should be energetically close, and they are on the MP2lMP4 surface. All three DFT methods agree in placing structure 3c second in stability to 3a, but only the local density S-VWN method also places structure 3b close to those two. The S-VW method also agrees with the MP2 results in finding structure 3f to be least stable. 3f has three Na-0 bonds and no hydrogen bonds. Structures 3d and 3e lie intermediate in energy on the MP2 surface between the nearly degenerate 3b and 3c and structure 3f. Isomer 3d has three Na-0 bonds and one hydrogen bond while 3e is another cyclic structure, similar to 3a, but with two Na--0 and two hydrogen bonds. None of the DFT methods finds structures 3d and 3e to be especially close in energy, S-VWN placing them closest. Overall, the S-VWN method orders the isomers in closest agreement with MP2, while the actual magnitudes of the energy separations are better replicated by B3-LYP and B-LYP. Because of the complexity of the Na(HZO)3 surface, a comprehensive search for false minima, such as was undertaken for Na(H?O)?, was not feasible. Such a study is, however, significant in assessing the suitability of using DFT techniques to describe hydrogen bonded systems, and is being considered.

4. Conclusions The results presented here for the Na(H20)n (n = l-3) structures indicate that the three DFT-based methods predict molecular geometries well, even for compounds whose structures depend on hydrogen

94

D. Bacelo. Y. IshikawalJoumal

of Molecular Shucture

bonding. All three DFT methods proved able to accurately locate global minima on the convoluted potential surfaces of these hydrogen bonded systems. Overall, the B3-LYP potential produced geometries which more consistently agree with the MP2 geometries than did either of the other two potentials, but even the local density S-VWN gave credible results. However, the energy orderings of the local minima of Na(H20)3 predicted by all the DFT methods differ from the MP2 and MP4 order. Moreover, even the DFT methods which include nonlocal exchange corrections failed to place isomer pairs which are structurally similar similarly close in energy. In this regard the LDA method proved superior. In addition, two of the DFT methods, B-LYP and S-VWN, displayed a tendency to introduce spurious minima into the Na(H20)* potential surface. The solution to the drawbacks of the DFT techniques awaits the development of systematic schemes for the incremental improvement of the exchange-correlation functional. The development of such schemes depends, in turn, upon a more complete understanding of the physics of the density functional formalism. Until such advances in understanding take place, it will remain possible to incorporate into DFT necessary but not sufficient conditions for accurate description of molecular electronic structures.

Acknowledgements This work has been supported by the NSF through the EPSCoR program (Cooperative Agreement No. OSR-9452893). The authors gratefully acknowledge the assistance of Dr R.C. Binning, Jr., of the NIH RCMI Center for Molecular Modeling and Computational Chemistry at the University of Puerto Rico in preparing this manuscript.

References [1] (a) P. Jena, S.N. Khanna,

B.K. Rao (Eds.), Physics and Chemistry of Finite Systems: from Clusters to Crystals, ~01s. 1 and 2, NATO AS1 Series, v. C374, Kluwer, Dordrecht, 1992. (b) H. Nakatsuji, H. Nakai, M. Hada, in: D.R. Salahub, N.

(Theochem) 425 (1998) 87-94

Russo (Eds.), Metal-Ligand Interactions: from Atoms, to Clusters, to Surfaces, Kluwer, Dordrecht, 1992, pp. 25 l-285. [2] C.P. Schulz, R. Haugstltter, H.-U. Tittes, I.V. Hertel, 2. Physik D 10 (1988) 279. [3] I.V. Hertel, C. Hiiglin, C. Nitsch, C.P. Schulz, Phys. Rev. Lett. 67 (1991) 1767. [4] K. Fuke, K. Tsukamoto, M. Sanekata, in: P. Jena, S.N. Khanna, B.K. Rao (Eds.), Physics and Chemistry of Finite Systems: from Clusters to Crystals, vol. 2, NATO AS1 Series, v. C374, Kluwer, Dordrecht, 1992, p. 925. [5] F. Misaizu, K. Tsukamoto, M. Sanekata, K. Fuke, Chem. Phys. Len. 188 (1992) 241. [6] L.A. Curtiss, E. Kraka, J. Gauss, D. Cremer, J. Phys. Chem. 91 (1987) 1080. 17) K. Hashimoto, S. He, K. Morokuma, Chem. Phys. Lett. 206 (1993) 297. [S] R.N. Bamett and U. Landman, Phys. Rev. Lett. 70 (1993) 1775. [9] Y. Ishikawa, R.C. Binning, Jr., H. Sekino, Int. J. Quantum Chem. S29 (1995) 669. [lo] C. Moller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [1 1] J.A. Pople, H.B. Schlegel, R. Krishnan, D.J. DeFrees, J.S. Binkley, M.J. Frisch, R.A. Whitesides, R.J. Hout, W.J. Hehre, Int. J. Quantum Chem. Symp. S 15 (198 I) 269. (121 L.A. Curtiss, K. Raghavachari, in: S. Langhoff (Ed.), Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy, Kluwer, Dordrecht, 1994. [13] (a) T. Ziegler, Chem. Rev. 91 (1991) 651.(b) P. Mlynarski, D.R. Salahub, J. Chem. Phys. 95 (1991) 6050. [14] B.G. Johnson, P.M.W. Gill, J.A. Pople, J. Chem. Phys. 98 (1993) 5612. [15] J.P. Perdew, Y. Wang, Phys. Rev. B 45 (1992) 13244. [16] AD. Becke, J. Chem. Phys. 88 (1988) 1053. [17] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [ 181 M. Kieninger, S. Suhai, Int. J. Quantum Chem. 52 (1994) 465. [19] J.J. Novoa, C. Sosa, J. Phys. Chem. 99 (1995) 15837. [20] (a) A.D. Becke, J. Chem. Phys. 98 (1993) 1372. (b) A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [21] P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, J. Phys. Chem. 98 (1994) 1 I 623. [22] 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. AlLaham, V.G. Zakrzewski, J.V. Ortiz, J.B. Foresman, J. Cioslowski, B.B. Stefanov, A. Nanayakkara, M. Challacombe, 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. Head-Gordon, C. Gonzalez, J.A. Pople, GAUSSIAN 94, Gaussian, Inc., Pittsburgh, PA, 1995. [23] W.J. Hehre, R. Ditchtield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. [24] E. Honegger, S. Leutwyler, J. Chem. Phys. 88 (1988) 2582. [25] J.C. Slater, Phys. Rev. 81 (1951) 385. [26] S.J. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200.