Consequences of accidental degeneracy within density functional theory: the enigmatic structure of boron nitrosyl

15 April 2002

Chemical Physics Letters 356 (2002) 7–13 www.elsevier.com/locate/cplett

Consequences of accidental degeneracy within density functional theory: the enigmatic structure of boron nitrosyl Wei-Quan Tian, Galina Orlova, John D. Goddard

*

Department of Chemistry and Biochemistry, University of Guelph, Guelph, Ont., Canada N1G 2W1 Received 14 November 2001; in final form 4 February 2002

Abstract Boron nitrosyl was examined with DFT, MP2, CCSD(T), and CASSCF. Poor overlap of the atomic orbitals on boron with the molecular orbitals of NO causes an accidental near degeneracy of the highest occupied r- and porbitals. A near degeneracy of three electronic states results. DFT locates all three electronic states but finds unstable solutions for the 3 R
and 3 P linear isomers. CASSCF/6-311+G(3df) and MP2/6-311+G(2df) predict 3 R
and 3 P isomers but do not locate the 3 A00 bent structure. CCSD(T)/6-311+G(d), predicts the 3 A00 bent isomer to be 1.7 and 4.2 kcal/mol lower than the 3 R
and 3 P linear isomers, respectively. Ó 2002 Elsevier Science B.V. All rights reserved.

1. Introduction The formation and reactivity of the XNO–NXO–XON molecules (X ¼ B, Al, Ga, In, Tl) have attracted the focused attention of a number of experimental and theoretical chemists [1–10] due to growing interest in the chemistry of nitric oxide [8], and especially in the design of new materials such as ultrathin gate dielectrics [9]. Gallium and indium containing catalysts are used in the reaction between nitrogen oxide and methane [10]. Nitrosyls of the heavier elements, Al, Ga, In, and Tl, have been trapped in solid

*

Corresponding author. E-mail address: [email protected] (J.D. Goddard).

argon matrices and studied by IR spectroscopy by the Andrews’ group [2,3]. DFT and MP2 theoretical methods have been used for the assignment of the IR spectra. These combined experimental and theoretical studies showed that the XNO structures are lower in energy than the NXO and NOX isomers. Nitrosyls of Al, Ga, and In have been predicted to have linear structures with a 3 R
ground electronic state. The structure of TlNO has been the subject of discussion [3] since the B3LYP functional predicted a linear structure while the MP2 method yielded a bent form. In a recent paper [4] it was shown that the higher-level coupled-cluster method, CCSD(T), predicts a linear structure for TlNO similar to the GGA and hybrid DFT approaches.

0009-2614/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 ( 0 2 ) 0 0 2 5 3 - 1

8

W.-Q. Tian et al. / Chemical Physics Letters 356 (2002) 7–13

For the BNO–NBO–BON system, no experimental results have been reported as yet. According to previous theoretical investigations [5–7], boron nitrosyl is unique among the XNO (X ¼ B, Al, Ga, In, Tl) series of the group 13 elements in that BNO is a high-energy isomer, in contrast to the nitrosyls of the heavier elements. The most recent DFT and MP2 studies by Zhang and Zhou [7] have shown that the BNO isomer forms in an initial step on the B + NO reaction pathway. Subsequent isomerization yields the insertion product, NBO, which is ca. 70–80 kcal/mol lower in energy than BNO. Boron isonitrosyl, BON, is an even higher-energy isomer (ca. 90 kcal/mol above NBO). Triplet ground electronic states were assigned to all three isomers. The structure of BNO was not predicted unambiguously. The MP2/ 6-311+G(d) method yielded a 3 A00 bent isomer. However, this bent form converged to a 3 R
linear structure during geometry optimization when a larger basis set, 6-311+G(2df), with f-type polarization functions, was employed. In contrast, GGA and hybrid DFT with both these basis sets predicted a 3 A00 bent isomer while the 3 R
linear form is a transition structure. It was suggested that higher-level theoretical methods should be employed. Two resonance structures which have been proposed [6] for the linear and bent forms of BNO are shown in Scheme 1. The DFT and MP2 geometries reported by Zhang and Zhou [7] resemble structures 1a and 2a, with an essentially single BN bond. We performed theoretical studies on boron nitrosyl using DFT, MP2, CCSD(T), and CASSCF methods to predict the ground electronic state as well as to examine the source of the different behaviors of these theoretical methods in structural predictions on BNO.

2. Computational details All predictions were made with GA U S S I A N 98 [11]. Hybrid DFT which incorporates delocalized Hartree–Fock exchange and pure local DFT have a different balance of dynamical and nondynamical electron correlation effects and, therefore, often perform differently [4]. Thus, both these DFT ap-

Scheme 1.

proaches were employed, within the unrestricted approximation. For pure local DFT within the generalized gradient approximation (GGA), the local B exchange functional [12] with the LYP [13] and P86 [14] correlation functionals (BLYP and BP86) were employed. The hybrid exchange functional B3 [15] with the correlation functional LYP (B3LYP) was used for the hybrid DFT predictions. The pruned (99,590) grid was used for numerical integration with DFT. Møller–Plesset theory to second order, MP2 [16] and the coupledcluster method, CCSD(T) [17], also were employed. All the post-Hartree–Fock predictions were performed with frozen core, within the unrestricted approximation. For a preliminary complete active space self-consistent field CASSCF [18–23] study, the active space included eight electrons and eight orbitals, five occupied and three virtual. This 8o/8e active space spans the entire p-system and the highest occupied and weakly bonding r-orbital. The lower strongly bonding r-orbitals predominantly constructed from the 2s and 2pz atomic orbitals were excluded from active space. Triply split valence basis sets augmented with diffuse functions and with d-type and f-type polarization functions, 6-311+G(d), 6-311+G(2df) and 6-311+G(3df) [24,25], were employed. The stabilities of the unrestricted DFT solutions were checked using the tests on stability [26] implemented in GA U S S I A N 98. Harmonic vibrational frequencies and infrared intensities at the DFT, MP2, and CASSCF levels were determined via analytic second derivative techniques [27,28]. For the CCSD(T) method, gradients and harmonic vibrational frequencies were determined by first and second numerical differences of the energies.

W.-Q. Tian et al. / Chemical Physics Letters 356 (2002) 7–13

3. Results and discussion 3.1. DFT with the hybrid B3LYP and GGA BLYP and BP86 exchange-correlation functionals A 3 R linear structure of BNO and the lower energy 3 A00 bent isomer, which resemble structures 1a and 2a, respectively, were found. The C1v 3 R
structure is a saddle point with the hybrid B3LYP functional in agreement with earlier B3LYP/6311+G(d) predictions [6,7]. The two degenerate imaginary frequencies correspond to bending modes. In contrast to previously reported results [6,7], both GGA functionals predict the C1v 3 R
structure to be a minimum rather than a transition structure. Another linear isomer, which resembles structure 1b, with a BN double bond, was located with both the GGA and hybrid DFT approaches. This isomer possesses a 3 P electronic state and is slightly lower in energy than the 3 R
linear structure. For linear BNO, the KS solutions have a tendency to converge to a 3 R
electronic state.

9

Thus, to obtain the 3 P isomer, the starting geometry should be considerably changed from that of the optimized 3 R
structure, to a short B@N
and a long N@O bond of ca. bond of ca. 1.19 A
1.4 A in the initial geometry for the 3 P electronic state. Geometries, energetics, and harmonic vibrational spectra predicted for the 3 R
; 3 P; and 3 00 A isomers with the hybrid DFT and GGA functionals are presented in Table 1. Analysis of the KS wave functions shows that the 3 R
and 3 P electronic states predicted for linear BNO with the DFT functionals in fact are not the true ground electronic states. The KS orbitals for the 3 R
; 3 P, and 3 A00 electronic states are illustrated in Fig. 1. Both the GGA BP86 and BLYP functionals and the hybrid B3LYP functional yield an unconventional ordering of the pand r-orbitals for the 3 R
ground electronic state, i.e., the ‘unpaired’ electrons are in the p-orbitals which are not the highest energy occupied orbitals. the 3 P electronic state also has an unusual ordering of orbitals, i.e., the formally singly occupied

Table 1
and degrees; total energies, Et , in hartree; relative energies, DEt , in kcal/mol (values in parentheses include ZPVE Geometries, in A corrections); harmonic vibrational frequencies, in cm
1 ; along with IR intensities, in km/mol (values in parentheses); predicted for the 3
3 R ; P, and 3 A00 structures of BNO with the B3LYP, BLYP, and BP86 functionals using the 6-311+G(2df) basis set Method/system

Geometries

Et

DEt

Frequencies

B3LYP 3
R

rB–N ¼ 1:412; rN–O ¼ 1:184

)154.695284 ()154.688945) )154.697198 ()154.688573) )154.699718 ()154.693042)

2.8 (2.6) 1.6 (2.8) 0.0 (0.0)

p ¼ 632i ð312  2Þ; r ¼ 931 ð225Þ; r ¼ 1852 ð1Þ p ¼ 254 ð43Þ; p ¼ 415 ð2Þ; r ¼ 1207 ð83Þ; r ¼ 1911 ð27Þ a0 ¼ 327 ð17Þ; a0 ¼ 925 ð22Þ; a0 ¼ 1678 ð332Þ

)154.678635 ()154.668415) )154.681332 ()154.672285) )154.687032 ()154.680718)

5.3 (7.7) 3.6 (5.3) 0.0 (0.0)

r ¼ 871 ð203Þ; p ¼ 949 ð745  2Þ; r ¼ 1717 ð8Þ p ¼ 401 ð2Þ; p ¼ 574 ð185Þ; r ¼ 1158 ð78Þ; r ¼ 1839 ð49Þ a0 ¼ 331 ð5Þ; a0 ¼ 913 ð46Þ; a0 ¼ 1528 ð465Þ

)154.696172 ()154.686199) )154.701676 (154.692754) )154.705181 ()154.698776)

5.6 (7.9) 2.2 (3.8) 0.0 (0.0)

r ¼ 871 ð204Þ; p ¼ 877 ð495  2Þ; r ¼ 1752 ð1Þ p ¼ 416 ð2Þ; p ¼ 454 ð32Þ; r ¼ 1183 ð77Þ; r ¼ 1864 ð68Þ a0 ¼ 333 ð5Þ; a0 ¼ 907 ð88Þ; a0 ¼ 1572 ð575Þ

3

P

rB–N ¼ 1:266; rN–O ¼ 1:216

3

A00

rB–N ¼ 1:394; rN–O ¼ 1:196 \BNO ¼ 151:5

BLYP 3
R

rB–N ¼ 1:432; rN–O ¼ 1:201

3

P

rB–N ¼ 1:282; rN–O ¼ 1:228

3

A00

rB–N ¼ 1:382; rN–O ¼ 1:217 \BNO ¼ 152:5

BP86 3
R

rB–N ¼ 1:436; rN–O ¼ 1:194

3

P

rB–N ¼ 1:285; rN–O ¼ 1:219

3

A00

rB–N ¼ 1:377; rN–O ¼ 1:209 \BNO ¼ 154:2

10

W.-Q. Tian et al. / Chemical Physics Letters 356 (2002) 7–13

Fig. 1. The predicted order of the KS orbitals for the 3 R
; 3 P, and 3 A00 electronic states of BNO and for the 3 R
electronic state of AlNO with DFT.

r-orbital is not the highest energy orbital. Such an ‘incorrect’ ordering of orbitals may cause an instability of the KS solutions. Indeed, tests on stability show that all the 3 R
and 3 P solutions, except for the 3 R
solution from B3LYP, are unstable. The re-optimization of the SP KS solutions in the space of the occupied orbitals without any symmetry constraints yields a small (ca. 0.1 kcal/mol) lowering in the energy. The results of

stability tests are reported in Table 2. The stability vectors correspond to the rotation of the p- and rorbitals and, for the 3 P solutions, lead to the ground 3
R electronic state with a conventional ordering of orbitals. In turn, for the 3 R
solutions, stability vectors lead to the ground 3 P electronic state with a conventional ordering of orbitals. However, we were not able to locate these states with DFT. The 3 A00 bent isomer was found to be the lowest in energy with DFT. Predicted geometries correspond more closely to the resonance structure, 2a, with a B–N single bond. The formally singly occupied p-orbitals are the highest energy orbitals as is typical of a ground 3 A00 electronic state and the SP KS solutions are stable. Interestingly, the 3 R
solution from the B3LYP functional is stable however the 3 R
C1v structure corresponds not to a minimum but to a transition structure. Such an alternative, i.e., a stable KS solution for a saddle point structure in contrast to a minimum with an unstable KS solution has been þ observed previously for the 2 Rþ u C5 radical cation. For this system, DFT predicts an ‘incorrect’ ordering of the KS orbitals, similar to the 3 R
BNO case. The BP86 functional generates a stable KS þ solution and predicts 2 Rþ u C5 to be a transition structure while the hybrid B3LYP functional generates an unstable 2 Rþ u solution which corresponds to a minimum [29]. The analysis of the KS orbitals reveals that the near degeneracy of the 3 R
and 3 P states occurs due to an accidental near-degeneracy of the p-orbitals and the highest occupied r-orbital. The

Table 2 Eigenvalues (in hartrees) and eigenvectors of the stability matrices for the 3 R
and 3 P unstable solutions generated for the C1v structures of BNO with the B3LYP, BLYP, and BP86 functionals using the 6-311+G(2df) basis set System/method

Root

Eigenvalue

3

P/B3LYP R
=BLYP

1 1, 2

)0.0000158 )0.0043263

3

P=BLYP

1

)0.0044093

3

R
=BP86

1, 2

)0.0055296

3

3

P=BP86

1

)0.0000296

Eigenvector 0.992 0.259 0.954 0.117 0.984 0.195 )0.190 0.680 )0.661 0.997

The numbering and order of the highest orbitals are shown in Fig. 1.

Orbital rotation

El. st.

b : 9p ! 11r a : 11r ! 13p b : 9r ! 10p a : 9r ! 13p b : 9p ! 11r a : 11r ! 12p a : 11r ! 13p b : 9r ! 10p b : 9r ! 11p b : 9p ! 11r

3

R

3

P

3

R

3

P

3

R

W.-Q. Tian et al. / Chemical Physics Letters 356 (2002) 7–13

11

and 3 A00 structures of BNO found by DFT. The 8o/ 8e active space, based on natural orbitals, spans the ð7r2 1p2 2p2 3p1 4p1 5p0 6p0 8r0 Þ orbitals of the 3 R
electronic state, the ð6r2 2p2 3p2 7r1 4p1 5p0 6p0 8r0 Þ orbitals of the 3 P electronic state, 2 2 2 1 1 0 0 0 and the ð7a0 1a00 8a0 9a0 2a00 10a0 3a00 11a0 Þ 00 orbitals of the 3 A electronic state. Therefore, the p-system and the highest occupied weakly bonding r-orbital which causes the degeneracy problems are included. The predicted results are reported in Table 3. Only two linear isomers corresponding to the 3 P and 3 R
electronic states were located with CASSCF. These two linear isomers are nearly isoenergetic, with an energy difference of only 0.02 kcal/mol in favor of the 3 P isomer. The inclusion of ZPVE corrections causes the 3 R
isomer to be 0.77 kcal/mol lower in energy than the 3 P isomer. The 3 A00 bent isomer was not located since the starting geometry corresponding to the bent structure converged to the 3 P linear isomer during geometry optimization. The analysis of the electron configurations shows that for both isomers, the ground state electronic configuration is dominant with a weight of ca. 90%. The lack of multireference character indicates that for BNO a single reference based method such as CCSD(T) should predict accurate results.

poor overlap of the 2s atomic orbital on boron with the r-system of NO accounts for this unusual p–r degeneracy. The BP86/6-31G(d) predictions of the chemically closest analogue, AlNO, confirm this supposition. For AlNO, the highest occupied r-orbital has well-developed Al–N and N–O bonding character and is notably lower in energy than the p-orbitals (see Fig. 1). As a result, the 3 R
electronic state is a true ground state, with a conventional ordering of orbitals and the 3 R
KS solution is stable. For boron nitrosyl, it is impossible to avoid the p–r degeneracy within C1v symmetry and, thus, the hybrid B3LYP and GGA BLYP and BP86 functionals are bound to generate unstable solutions. The stability tests yield only very small lowering in energy with the re-optimization of the SP KS solutions in the space of the occupied orbitals without any symmetry constraints. Stable solutions, if found, could not give notable differences in energetics and geometries compared to those found with the unstable solutions. However, the relative energetics may change due to the very delicate energy differences between the 3 R
; 3 P, and 3 A00 electronic states. 3.2. CASSCF(8, 8) The CASSCF approximation is one remedy to ‘sort out’ nearly degenerate electronic states. However, as it largely neglects dynamical electron correlation, CASSCF might not be the best method in the case of BNO. Boron nitrosyl appears not to be a ‘multireference’ system since the near degeneracy occurs accidentally due to a poor interaction between the boron atom and the NO fragment rather than to any spatial symmetry problems. We performed a preliminary CASSCF(8, 8)/6-311+G(3df) study of the 3 R
; 3 P,

3.3. MP2 and CCSD(T) The unrestricted Hartree–Fock (UHF) wave functions which were used as reference functions for the MP2 and CCSD(T) show the same types of instabilities for the 3 R
and 3 P electronic states as do the KS solutions. The instabilities in UHF which do not involve spatial symmetry breaking are not expected to affect MP2 and

Table 3
and degrees; total energies, Et , in hartree; relative energies, DEt , in kcal/mol (values in parentheses include ZPVE Geometries, in A correction); and harmonic vibrational frequencies, in cm
1 predicted for BNO isomers with the CASSCF(8, 8)/6-311+G(3df) method El. st.


Geometries

Et

DEt

Frequencies

rB–N ¼ 1:387; rN–O ¼ 1:190

)153.989006 ()153.981782) )153.989042 ()153.980547)

0.02 ()0.77) 0.0 (0.0)

p ¼ 164; r ¼ 995; r ¼ 1849 p ¼ 300; p ¼ 417; r ¼ 1135; r ¼ 1876

3

R

3

P

rB–N ¼ 1:270; rN–O ¼ 1:189

3

A00

Geometry optimization converges to 3 P

12

W.-Q. Tian et al. / Chemical Physics Letters 356 (2002) 7–13

CCSD(T) predictions. The results predicted for BNO with the MP2 and CCSD(T) methods are reported in Table 4. The MP2/6-311+G(d) method yields nearly isoenergetic 3 R
; 3 P, and 3 00 A isomers, with an insignificant energetic preference for the 3 R
structure. The inclusion of the ZPVE correction does not change qualitatively the relative stabilities. The geometries of the 3 R
; 3 P; and 3 A00 isomers correspond to the structures 1a, 1b, and 2a, respectively. An inclusion of f-type polarization functions in the basis set does qualitatively change the MP2 predictions. An initial bent geometry leads to a linear form during geometry optimization in accord with earlier results with MP2 including f-functions [7]. However, this effectively linear form, corresponds to the 3 P isomer, with a BN double bond rather than to the 3 R
isomer. With the 6-311+G(2df) basis set, the 3 R
isomer remains lower in energy than the 3 P one. The higher-level treatment of electron correlation by CCSD(T) frequently predicts very accu-

rate results for the systems which do not possess significant multireference character. The CCSD(T)/6-311+G(d) method yields the 3 R
1a and 3 A00 2a isomers, with a BN single bond and the 3 P 1b isomer, with a BN double bond. The bent isomer is the lowest in energy, while the 3 P isomer is the highest one. Unlike MP2, an inclusion of f-type polarization functions does not affect qualitatively the predictions from CCSD(T). The CCSD(T)/6-311+G(2df) method predicts the 3 A00 bent isomer with a BN single bond to be a minimum on the potential energy surface. CCSD(T) overcomes the instability from which the DFT functionals and UHF reference functions suffer. However, for the 3 R
isomer, there is a variational collapse into the low-lying 3 A00 electronic states which results in harmonic vibrational frequencies for the bending which are too large and to too large a ZPVE correction. Therefore, ZPVE corrections should not be included in the relative energies from CCSD(T).

Table 4
and degrees; total energies, Et , in hartree; relative energies, DEt , in kcal/mol (values in parentheses include ZPVE Geometries, in A correction); and selected harmonic vibrational frequencies, in cm
1 ; along with IR intensities, in km/mol (values in parentheses), predicted for BNO isomers with the MP2 and CCSD(T) methods using the 6-311+G(d) and 6-311+G(2df) basis sets Method

El. st.

Geometries

Et

DEt

Frequencies

)2.4 ()0.5) 1.7 (3.3) 0.0 (0.0)

p ¼ 872 ð778  2Þ; r ¼ 951 ð208Þ; r ¼ 1813 ð8Þ p ¼ 444 ð3Þ; p ¼ 452 ð2Þ; r ¼ 1269 ð140Þ; r ¼ 2081 ð170Þ a0 ¼ 185 ð45Þ; a0 ¼ 960 ð171Þ; a0 ¼ 2004 ð9Þ

a

MP2/+G(d)

3

R


rB–N ¼ 1:416; rN–O ¼ 1:194

3

P

rB–N ¼ 1:257; rN–O ¼ 1:209

3

A00

rB–N ¼ 1:421; rN–O ¼ 1:195; \BNO ¼ 153:9

)154.271197 ()154.260925) )154.264603 ()154.254932) )154.267333 ()154.260165)

3

R
P

rB–N ¼ 1:411; rN–O ¼ 1:190 rB–N ¼ 1:253; rN–O ¼ 1:204

)154.345488 )154.340765

0.0 3.0

3

R


rB–N ¼ 1:431; rN–O ¼ 1:194

3

P

rB–N ¼ 1:281; rN–O ¼ 1:222

3

A00

rB–N ¼ 1:434; rN–O ¼ 1:204 \BNO ¼ 148:3

)154.311286 ()154.297854) )154.307335 ()154.298813) )154.313944 ()154.307290)

1.7 (5.9) 4.2 (5.3) 0.0 (0.0)

rB–N ¼ 1:422; rN–O ¼ 1:198 \BNO ¼ 152:0

)154.389402 ()154.382711)

MP2/+G(2df) 3

CCSD(T)/+G(d)

CCSD(T)/+G(2df) 3 00 A a

For the basis set the common ‘6-311’ is omitted.

p ¼ 1086; r ¼ 1554; r ¼ 2170 p ¼ 274; p ¼ 398; r ¼ 1196; r ¼ 1872 a0 ¼ 281; a0 ¼ 927; a0 ¼ 1714 a0 ¼ 272; a0 ¼ 931; a0 ¼ 1733

W.-Q. Tian et al. / Chemical Physics Letters 356 (2002) 7–13

4. Conclusions For linear BNO, the highest occupied 2p- and 7r-orbitals are accidentally nearly degenerate due to the poor overlap of the 2s atomic orbital on boron with the r-orbitals of NO. As a result, the 3 R
; 3 P, and 3 A00 electronics states are nearly degenerate. DFT locates all the three electronic states but, being an ‘SCF-type’ method, generates unstable KS solutions for the 3 R
and 3 P linear structures. Thus, the caveat ‘Check the stability of your SCF solutions’ should be extended to DFT. The higherlevel CCSD(T) method also locates the 3 R
; 3 P, and 3 A00 structures and predicts the 3 A00 isomer to be lowest in energy in qualitative accord with DFT. The 3 P isomer is higher in energy than the 3 R
one with CCSD(T), in contrast to the DFT predictions. A more accurate evaluation of the relative energies for these linear isomers is hindered with DFT due to the instability of the KS solutions and with CCSD(T) due to variational collapse of the 3 R
structure. The CASSCF method which deemphasizes dynamical electron correlation and MP2 fail to locate the 3 A00 bent structure.

Acknowledgements Financial support of this research by the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged.

References [1] G.K. Ruschel, D.W. Ball, High Temp. Mater. Sci. 37 (1997) 63. [2] L. Andrews, M. Zhou, W.D. Bare, J. Phys. Chem. A 102 (1998) 5019.

13

[3] L. Andrews, M. Zhou, X. Wang, J. Phys. Chem. A 104 (2000) 8475. [4] G. Orlova, J.D. Goddard, Mol. Phys. 100 (2002) 483. [5] J.A. Harrison, R.G.A.R. Maclagan, Chem. Phys. Lett. 155 (1989) 419. [6] S. Su, THEOCHEM 430 (1998) 137. [7] L. Zhang, M. Zhou, Chem. Phys. 256 (2000) 185. [8] M.F. Perutz, Nature 380 (1996) 205. [9] E.P. Gusev, H.C. Lu, E. Garfunkel, T. Gustafsson, M.L. Green, D. Brasen, W.N. Lennard, J. Appl. Phys. 84 (1998) 2980. [10] E. Kikuchi, M. Ogura, I. Terasake, Y. Goto, J. Catal. 161 (1996) 465. [11] M.J. Frisch, et al., GA U S S I A N 98, Gaussian Inc., Pittsburgh, PA, 1998. [12] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [13] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [14] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [15] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [16] C. Møller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [17] K. Raghavachari, J.A. Pople, Int. J. Quantum Chem. 20 (1981) 167. [18] D. Hegarty, M.A. Robb, Mol. Phys. 38 (1979) 1795. [19] R.H.E. Eade, M.A. Robb, Chem. Phys. Lett. 83 (1981) 362. [20] H.B. Schlegel, M.A. Robb, Chem. Phys. Lett. 93 (1982) 43. [21] F. Bernardi, A. Bottini, J.J.W. McDougall, M.A. Robb, H.B. Schlegel, Faraday Symp. Chem. Soc. 19 (1984) 137. [22] N. Yamamoto, T. Vreven, M.A. Robb, M.J. Frisch, H.B. Schlegel, Chem. Phys. Lett. 250 (1996) 373. [23] M.J. Frisch, I.N. Ragazos, M.A. Robb, H.B. Schlegel, Chem. Phys. Lett. 198 (1992) 524. [24] A.D. McLean, G.S. Chandler, J. Chem. Phys. 72 (1980) 5639. [25] R. Krishnan, J.S. Binkley, R. Seeger, J.A. Pople, J. Chem. Phys. 72 (1980) 650. [26] R. Bauernschmitt, R. Ahlrichs, J. Chem. Phys. 104 (1996) 9047. [27] B.G. Johnson, M.J. Frisch, Chem. Phys. Lett. 216 (1993) 133. [28] B.G. Johnson, M.J. Frisch, J. Chem. Phys. 100 (1994) 7429. [29] G. Orlova, J.D. Goddard, J. Phys. Chem. A (submitted).