- Email: [email protected]
An ab initio and DFT study of (N2)2 dimers
4 June 1999
Chemical Physics Letters 306 Ž1999. 71–77
An ab initio and DFT study of žN2 /2 dimers O. Couronne ) , Y. Ellinger Laboratoire d’Etude Theorique des Milieux Extremes, Ecole Normale Superieure, 24 rue Lhomond, F-75231 Paris Cedex 05, France ´ ˆ ´ Received 7 December 1998; in final form 6 April 1999
Abstract The structure of van der Waals dimers ŽN2 . 2 is studied using ab initio and density functional calculations. The potential energy surfaces corrected a priori for basis set superposition errors are necessary for determining both the geometries and the vibrational frequencies. With both ab initio MP2, MP4 and DFT PW91–PW91 levels of theory, the T-shaped and canted conformations appear to be the most stable, within 1–5 cmy1 of each other. The DFT PW91–PW91 dissociation energy is 67 cmy1 and an upper limit to the barrier to internal motion is 30 cmy1, both in excellent agreement with the values deduced from IR measurements. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction A large part of the mass of the universe is concentrated in dust particles. Due to the low temperature Ž10–50 K. of a number of celestial objects Žcold molecular clouds, comets, rings of giant planets, etc.. these particles behave as cold surfaces on which light molecules condense, forming what is known as ice mantles. This generic term, besides water, covers all types of organic or non-organic materials such as CH 4 , NH 3 , CH 3 OH, CO, CO 2 , or N2 whose accretion at surfaces is at the origin of the growth of interstellar grains. Before dealing with the structure of molecular solids and the possible chemistry that takes place at their surface under the extreme conditions prevailing in space, we found it necessary to consider the primary process of accretion, namely the formation of isolated dimers. )
Corresponding author. Present address: Department of Physical Chemistry, University of Geneva, 30 quai Ansermet, CH-1211 Geneva 4, Switzerland. E-mail: [email protected]
The first dimer, the study of which is reported in this Letter, is ŽN2 . 2 ; it may play an important role in the atmospheric chemistry of Neptune and Titan. Experimental evidence is limited but suggests that the most stable dimer has either a T-shaped or X-shaped geometry, with the centers of mass of the ˚ monomers separated by approximately 3.7–4.2 A w1,2x. From a theoretical point of view, the ŽN2 . 2 dimers have been studied with ab initio w3–8x and model potential w9–12x calculations. Most of these studies have predicted up to four local minima ŽTshaped, X-shaped, parallel and canted. but the energetic ordering varies with the method. The Møller–Plesset level of theory is usually accepted as a good starting point to account for the correlation effects required to describe weak interactions in van der Waals systems w13–15x. We will use the results obtained with this methodology for comparison with DFT calculations on this type of compounds. Although DFT gives remarkable results for solid structures as well as for covalent and ionic systems, it is not completely clear yet if it can
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 4 3 1 - 5
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
72
systematically deal with van der Waals supermolecules w16,17x even if clear-cut results have been obtained for rare gas dimers w18–20x, complexes between molecular nitrogen and benzene w21x and Cl 2 –C 2 H 4 w20x. Beyond the test case of the ŽN2 . 2 dimer, the present work is intended to probe the quality of DFT calculations in view of future studies on N2 , CO and CO 2 ices.
energy of the dimer in the molecular basis BD at geometry r D . The second term is a summation running over the energies of the monomers, each being calculated in the basis BD at their geometries in the dimer. The last term is a summation running over the energies of the monomers, calculated in the basis BM of each monomer at their geometries in the dimer. The BSSE-corrected PES is used to determine the geometry as well as the vibrational frequencies. The binding energy is given by:
2. Computational background D E s Ecorr y The ab initio calculations were performed at the MP2 and MP4 levels of theory as implemented in the GAUSSIAN94 package w22x. The standard 6-31 q G ) basis set was used. It has already proven reliable for these types of problems as well as for the description of molecular systems with a diffuse wavefunction such as negative ions w23x. Density functional theory ŽDFT. calculations were performed applying the Kohn–Sham ŽKS. formalism w24x. Two parametrizations of the exchange-correlation functionals which include gradient corrections within the generalized gradient approximation ŽGGA. were used: B3LYP w25,26x and PW91–PW91 w27x. DFT calculations were performed with GAUSSIAN98 w28x. In the DFT calculations, the atomic orbitals were constructed from the following contraction patterns: Ž7111r411r1) . w29x. This basis set is the closest to the basis used in ab initio calculations as shown by the comparison made with the results previously obtained on weakly interacting species w21x. The grids considered for the computations comprised 128 radial shells. At the end of the SCF procedure, each radial shell had 302 angular points. One of the main characteristics of the present study is that the considered potential energy surface ŽPES. is corrected for the basis set superposition error ŽBSSE. according to the counterpoise method w30x. The corrected energy is obtained as: Ecorr s ED Ž BD ,r D . y
Ý Ž EM Ž BD ,r D . Ms1,2
yEM Ž BM ,r D . . ,
Ž 1.
where the indices M and D refer to the monomer and dimer, respectively. The first term stands for the
Ý
EM Ž BM ,r M . .
Ž 2.
Ms1,2
It is now well established that this correction is essential for the determination of molecular structures bound by van der Waals interactions w31x. But unlike many studies where the BSSE correction is applied to the equilibrium structure, i.e. after the minimization has been completed, it is applied here to the entire optimization process. Therefore, we obtain a different surface, the minimum of which is the correct equilibrium geometry. We will show that, in these compounds where cohesion is dependent on weak intermolecular forces, neglecting BSSE or applying an inadequate correction can lead to dramatically erroneous results. As the BSSE-corrected energy is a sum of energies, it is possible to define a corrected Hessian H corr as the matrix of the second derivatives of the BSSEcorrected energy with respect to the nuclear coordinates, namely w32x: H corr s H D Ž BD ,r D . y
Ý ŽH M Ž BD ,r D . Ms1,2
yH M Ž BM ,r D . . ,
Ž 3.
where the terminology is the same as for the corrected energy and H indicates a Hessian matrix. Each second derivative matrix can be computed either analytically or numerically. Frequencies and normal modes corrected for BSSE are then obtained by diagonalization of the force constant matrix Fcorr which is derived from H corr . In such a way, the zero-point energy ŽZPE. obtained from these frequencies is also corrected for BSSE.
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
73
3. Results The geometry of the dimer in any conformation investigated can be defined by four parameters, namely the distance between the centers of mass R, two angles u 1 and u 2 , and one dihedral angle w ŽFig. 1.. Six types of conformations of ŽN2 . 2 have been studied: linear, parallel, T-shaped, X-shaped, canted and crystal ŽFig. 2., constrained to D`h , C2 Õ , D2 h , D2 d and C1 symmetries, respectively. 3.1. Symmetry constrained optimization The first step in our study consists in an optimization along the R coordinate. In these calculations the internuclear distance in the N2 molecule is the same for the dimer as for the monomer. The position of the minimum R 0 and the binding energy Emin are obtained by fitting n-order Ž n s 2–10. polynomials to 3–11 points. The BSSE correction appears to be essential for a reliable evaluation of both the equilibrium geometry and ZPE ŽFig. 4a.. More precisely, Table 1 shows that R min is systematically shifted toward longer distances, with a ˚ according to the correction varying from 0.2 to 0.4 A type of dimer. The smaller corrections appear in the DFT treatment. It can also be seen that the binding energy is systematically larger with the PW91–PW91 functional than in ab initio MP2 and MP4 calculations. In all cases it is greatly reduced when the BSSE correction is applied Žabout 80% for all conformations.. The BSSE-corrected potential energy curves obtained using MP2, MP4, DFT B3LYP and DFT PW91–PW91 are presented in Fig. 3 for the T-shaped conformation which will be shown below
Fig. 1. General representation for a dimer. M 1 and M 2 are the middle points of monomers 1 and 2. R is the distance between these two points. u 1 and u 2 are the angles between ŽM 1 M 2 . and the molecular axis. w is the dihedral angle between the two monomers.
Fig. 2. Different types of ŽN2 . 2 structures studied, as defined by the triplet Ž u 1 , u 2 , w ..
to be one of the two stable structures of equal energies. At the MP2 level, all dimers are found to be local minima, the T-shaped and canted conformations being the most stable ones with equivalent energies. Parallel and X-shaped conformations are found to be very close to each other, and the linear conformation Table 1 Equilibrium distances and binding energies along the R coordi˚. nate Ž D E in cmy1 and R min in A BSSE uncorrected
BSSE corrected
R min
DE
R min
DŽEqBSSE.
4.09 4.05 3.55 3.70 4.87
273 225 248 203 154
4.33 4.29 3.97 3.99 5.19
69 70 39 36 15
4.11 4.11 3.64 3.73 4.84
266 217 221 191 181
4.40 4.35 4.11 4.10 5.17
58 59 28 28 16
3.99 3.90 3.50 3.57 4.80
305 296 264 234 195
4.25 4.17 3.81 3.86 4.94
132 138 103 88 16
)
MP2r6-31qG T-shaped Canted X-shaped Parallel Linear MP4r6-31qG ) T-shaped Canted X-shaped Parallel Linear PW91–PW91r Ž7111r411r1) . T-shaped Canted X-shaped Parallel Linear
74
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
the dimer and will be discussed later when comparing larger N2 clusters to the crystal structure. 3.2. Minimum energy structures and Õibrational frequencies
Fig. 3. BSSE-corrected potential energy curves along the intermolecular coordinate R Žsee Fig. 1. for the T-shaped structure: MP2 Ž^., MP4 Žq., B3LYP Že. and PW91–PW91 Ž=..
is the least stable one. Binding energies at the MP4 level slightly differ from MP2, but do not interchange the relative stabilities of the dimers. For DFT calculations, all dimers were found to be non-bonded at the B3LYP level, implying that this functional is unable to reproduce dispersion forces for this type of molecular dimer. On the other hand, the PW91–PW91 results compare much better with the MP2 ones. The positions of the minima are close, ˚ toalthough systematically shifted by 0.10–0.15 A ˚ wards shorter distances Žexcept the 0.25 A shift for the linear dimer., the binding energy being a little larger. Very similar results have previously been obtained for the N2-benzene van der Waals complex w21x. The crystal conformation, as it appears in the a phase ŽFig. 2. is the third most stable arrangement after the T-shaped and canted conformations Ž D E s 56 cmy1 at the MP2 q BSSE level.. The distance between the two centers of mass optimized at the ˚ The MP2 q BSSE level of theory is R min s 4.2 A. experimental distance in the a phase of the crystal is ˚ Žcorresponding to a lattice parameter a s 5.7 4.0 A ˚A.. The a phase crystal conformation Ž u 1 s 90, u 2 s 35.3, w s 54.7. is not far from either the T-shaped Ž u 1 s 0, u 2 s 90, w s 0. or canted Ž u 1 s 45, u 2 s 45, w s 0. conformation. It can then be inferred that molecules tend to adopt T-shaped andror canted conformations but twist due to crystal constraints. This situation, however, has no physical meaning for
In order to assess the stability of the different dimers, frequencies have been computed on the corresponding MP2, MP4 and PW91–PW91 BSSE-corrected PES. As seen above, the force constant matrix corrected for BSSE is derived from a sum of Hessian matrices. Each Hessian can be computed either analytically or numerically on the MP2 q BSSE and PW91–PW91q BSSE surfaces, but only numerically for the MP4 q BSSE one. For the numerical calculations, we used a two-point differentiation. In order to obtain reliable results, especially for low frequencies, we had to use small displacements Ž d x . together with high accuracy in the energy value Ž d Eaccuracy .. We found, from comparison between analytical and numerical computations at the MP2 ˚ and level, that the parameters d x s 0.005 A d Eaccuracy s 10y7 a.u. give the best match between the two methods. MP4 q BSSE frequencies were then computed with these parameters. Frequencies have been obtained for the linear, parallel, X-shaped, canted and T-shaped conformations. All the frequencies of the T-shaped and the canted conformations are real, which confirms that these structures are two real minima ŽTable 2.. By Table 2 Internal modes Žcmy1 . of the potential surface for the T-shaped and the canted structures. Geometries are equilibrium geometries for each surface and force constant matrices were computed on the PESs corrected a priori for BSSE. The MP2 and PW91–PW91 force constant matrices were computed analytically, and the MP4 force constant matrix were computed numerically Žsee text for details.
T-shaped Bending Stretching Bending Out of plane Canted Bending Stretching Bending Out of plane
n MP 2
n MP4
n PW91 – PW91
b2 a1 b2 b1
;5 ;5 26 26
;5 ;5 16 21
26 27 37 39
b2 a1 b2 b1
;5 ;5 10 40
;5 ;5 8 32
22 24 27 61
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77 Table 3 Comparison between theoretical binding energies for the T-shaped and the canted conformations Ženergies in cmy1 .
T-shaped DŽEqBSSE. DŽEqBSSE.qZPE corr Canted DŽEqBSSE. DŽEqBSSE.qZPE corr Experiment Ž De . w1x
MP2
MP4
PW91–PW91
69 38
58 34
132 67
70 40
59 34
138 71 66
contrast, all the angular frequencies for the parallel, X-shaped, linear conformations are imaginary, which implies that these conformations are not stable.
75
Table 2 shows the vibrational frequencies for the internal modes of the most stable conformations: T-shaped and canted. ZPEs can be evaluated from these frequencies. The ZPE is a non-negligible quantity when compared to the binding energy: it represents up to 50%. The final binding energies are given in Table 3. Particularly good agreement is obtained between the DFT PW91–PW91 value and that of 66 cmy1 deduced from IR measurements by Long et al. w1x. 3.3. A priori Õersus a posteriori BSSE correction of binding energy One of the main points in this study is the a priori inclusion of the BSSE correction in the minimization process. This approach is seldom used, which is entirely justified when large energy differences are implied. In the present study, the energy differences involved are too small for this technique to be used. Using the T-shaped dimer as an example ŽFig. 4a. shows the result of an a posteriori BSSE correction on a minimum corresponding to an uncorrected PES. At the MP2r6-31q G ) level reported here, the error in the dissociation energy is 13 cmy1 , which represents 20% of the binding energy. The equilibrium distance between the centers of mass is under˚ In other cases, such as the estimated by 0.24 A.
Table 4 Comparison between internal modes Žcmy1 . before and after BSSE correction of the potential surface for the T-shaped structure ŽD 2 h .. Geometries are equilibrium geometries for each surface. The MP2 and PW91–PW91 force constant matrices are computed analytically, and the MP4 force constant matrix are computed numerically Žsee text for details.
Fig. 4. Stretching of the T-shaped conformation at the MP2r6-31 qG ) level. The dashed line represents the uncorrected energy and the plain line the energy corrected for BSSE. Ža. Difference between a posteriori correction of the BSSE Ž2. on the minimum of the non-corrected surface Ž1. and a minimization process on the BSSE-corrected surface Ž=.. Žb. ZPE corrections on the BSSEuncorrected and BSSE-corrected PES.
Before BSSE correction of the PES Bending Stretching Bending Out of plane After BSSE correction of the PES Bending Stretching Bending Out of plane
n MP 2
n MP4
n PW91yPW91
b2 a1 b2 b1
23 44 50 63
34 44 50 64
33 34 58 64
b2 a1 b2 b1
;5 ;5 26 26
;5 ;5 16 21
26 27 37 39
76
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
X-shaped conformation, the error in R increases up ˚ to 0.4 A. Using non-coherent ZPE is even worse. Fig. 4b shows two types of vibrational corrections, one computed on the corrected potential surface and hereafter referred to as ZPE corr Žtop., the other obtained from the uncorrected surface Žbottom.. Significant differences appear between the frequencies depending on whether they are corrected for BSSE or not ŽTable 4.. After correction, the low frequencies corresponding to the motion of one N2 molecule with respect to the other are lower than before correction. It was found that the stretching of the N2 bond is not affected by the formation of the dimer Žas anticipated from the conservation of the N–N bond length when going from the isolated monomer to the dimerized species.. It can be inferred from the curvature of the potential surfaces on Fig. 4a, that the second derivatives, and thus, the force constants, decrease when the BSSE correction is taken into account. By definition, BSSE is always negative. Since its magnitude decreases monotonically with intermolecular separation, it is quite obvious that correcting for BSSE results in a decrease of the PES curvature. A shallower surface implies lower frequencies for intermolecular modes and consequently a lower ZPE correction. For instance, at the MP2r6-31q G ) level, ZPE corr amounts to 31 cmy1 in the T-shaped dimer whereas the value directly obtained from the uncorrected PES is equal to 88 cmy1 which is 57 cmy1 higher. The present study clearly shows that erroneous results are obtained in such cases due to a critical overestimation of the vibrational frequencies. Table 5 Binding energies for the T-shaped conformation: comparison between a priori and a posteriori correction of the BSSE Ženergies in cmy1 ; =: instable.
DE a posteriori DEqBSSE DEqBSSEqZPE a priori DŽEqBSSE. DŽEqBSSE.qZPE corr
MP2
MP4
PW91– PW91
273
266
305
56 =
41 =
103 9
69 38
59 34
132 67
It can be seen in Table 5 that the T-shaped dimer is no longer bound at the MP2 and MP4 levels in these conditions Žthe same situation is found for the canted conformation.. 4. Concluding remarks In this contribution, we report a study of the structure and stability of the ŽN2 . 2 dimer at canted, T-shaped, X-shaped, parallel and linear structures plus that taken from the crystal. We have compared ab initio MP2 and MP4 results to the DFT ones. The present work provides further confirmation that the B3LYP functional is unable to properly describe the bonding forces in weakly interacting systems. By contrast, the exchange-correlation functional PW91– PW91 gives results which are not only close to the MP values but, even more satisfactorily, are in excellent agreement with IR measurements. The PW91–PW91 exchange-correlation functional, although still being semi-local and not correctly describing the Ry6 behaviour of the interaction energy at large distances, provides a very good approximation in the range of medium distances characteristic of the equilibrium structures of weakly interacting molecules. Moreover, this rather good performance of the PW91–PW91 functional, compared to that of B3LYP, can be attributed to its exchange part w21x. We have shown that PESs and ZPEs should be determined with explicit consideration of BSSE in the optimization procedure and calculations of the vibrational frequencies, a posteriori BSSE corrections possibly lead to erroneous equilibrium geometries, ZPE and binding energies. Under these conditions, the T-shaped and canted dimers are the only stable structures at the MP2, MP4 and DFT PW91– PW91 levels of theory, contrary to previous studies where equilibrium geometries were determined using PESs which were not BSSE corrected. It should be mentioned at this point that the same two conformations were found separated by 5 cmy1 in a recent study at the CCSDŽT. level with a posteriori BSSE corrections w33x. With proper a priori correction, the binding energy of the ŽN2 . 2 dimer is 38 cmy1 at the MP2 and MP4 levels of theory and 67 cmy1 for the DFT PW91–PW91 functional, in remarkable agreement with experiment.
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
The energy difference between the two stable conformations Ž1–6 cmy1 all corrections considered. is beyond the accuracy which can reasonably be expected when internal motions are treated at the harmonic level Žno coupling between vibrations.. No significant barrier could be found between the Tshaped and canted conformations, in agreement with the CCSDŽT. study w33x, which suggests that both conformations belong to the same potential well of the internal motions. The next saddle points, Xshaped and parallel can be seen as transition points for out-of-plane and in-plane Žnon-concerted. motions. We have verified that there is a direct path of decreasing energy connecting the parallel and Xshaped structures to the minima. These structures are about 30 cmy1 Ž43 K. higher than the equilibrium minima at all levels of theory. This value is in agreement with the 15–30 cmy1 Ž22–43 K. deduced from IR experiments w1x and is consistent with the permanent flip of the molecules in the crystal. Acknowledgements The authors would like to thank T.A. Wesołowski and Prof. N.C. Handy for stimulating remarks and suggestions in the course of this study. The financial support of CNRS Programme National de Planetolo´ gie ŽPNP. is greatly acknowledged. Part of the calculations were performed at IDRIS and CNUSC computing facilities. References w1x C.A. Long, G. Henderson, G.E. Ewing, Chem. Phys. 2 Ž1973. 485. w2x F. Carnovale, J.B. Peel, R.G. Rothwell, J. Chem. Phys. 88 Ž1988. 642. w3x F.H. Ree, N.W. Winter, J. Chem. Phys. 73 Ž1980. 322. w4x J. Tennyson, A. van der Avoird, J. Chem. Phys. 77 Ž1982. 5664. w5x P.J. Hay, R.T. Pack, R.L. Martin, J. Chem. Phys. 81 Ž1984. 1360. w6x H.-J. Bohm, R. Ahlrichs, Mol. Phys. 55 Ž1985. 1159. ¨ w7x F. Uhlik, Z. Slanina, A. Hinchliffe, J. Mol. Struct. ŽTHEOCHEM. 282 Ž1993. 271. w8x S. Raugei, G. Cardini, V. Schettino, Mol. Phys. 95 Ž1998. 477. w9x B. Jonsson, G. Karistrom, ¨ ¨ J. Chem. Phys. 74 Ž1981. 2896. w10x M.S.H. Ling, M. Rigby, Mol. Phys. 51 Ž1984. 855. w11x G. Brocks, A. van der Avoird, Mol. Phys. 55 Ž1985. 11.
77
w12x K.A. Franken, C.E. Dykstra, J. Phys. Chem. 97 Ž1993. 11408. w13x G. Granucci, M. Persico, Chem. Phys. Lett. 205 Ž1993. 331. w14x P. Hobza, H.L. Selze, E.W. Schlag, Chem. Rev. 94 Ž1994. 1767. w15x G. Chalasinski, M. Szczesniak, Chem. Rev. 94 Ž1994. 1723. ´ ´ w16x J.M. Perez-Jorda, ´ ´ A.D. Becke, Chem. Phys. Lett. 233 Ž1995. 134. w17x E. Ruiz, D.R. Salahub, A. Vela, J. Phys. Chem. 100 Ž1996. 12265. w18x Y. Zhang, W. Pan, W. Yang, J. Chem. Phys. 107 Ž1997. 7921. w19x D.C. Patton, D.V. Porezag, M.R. Pederson, Phys. Rev. B 55 Ž1997. 7454. w20x C. Adamo, V. Barone, J. Chem. Phys. 108 Ž1998. 664. w21x T.A. Wesołowski, O. Parisel, Y. Ellinger, J. Weber, J. Phys. Chem. A 101 Ž1997. 7818. w22x M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheeseman, T. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. Al-Laham, 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, Revision C.3, Gaussian Inc., Pittsburgh, PA, 1995. w23x O. Parisel, Y. Ellinger, C. Giessner-Prettre, Chem. Phys. Lett. 250 Ž1996. 178. w24x W. Kohn, L.J. Sham, Phys. Rev. A 140 Ž1965. 1133. w25x A.D. Becke, J. Chem. Phys. 98 Ž1993. 5648. w26x C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 Ž1988. 785. w27x J.P. Perdew, Y. Wang, in: P. Ziesche, H. Eschrig ŽEds.., Electronic Structure of Solids’ 91, Academie Verlag, Berlin, 1991. w28x M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, GAUSSIAN 98, Revision A.1, Gaussian Inc., Pittsburgh, PA, 1998. w29x F. Sim, D.R. Salahub, S. Chin, M.J. Dupuis, J. Chem. Phys. 95 Ž1990. 4317. w30x S.F. Boys, F. Bernardi, Mol. Phys. 19 Ž1970. 553. w31x P. Hobza, R. Zahradnık, ´ Chem. Rev. 88 Ž1988. 871. w32x S. Simon, M.D.J.J. Dannenberg, J. Chem. Phys. 105 Ž1996. 11024. w33x A. Wada, H. Kanamori, S. Iwata, J. Chem. Phys. 109 Ž1998. 9434.
Chemical Physics Letters 306 Ž1999. 71–77
An ab initio and DFT study of žN2 /2 dimers O. Couronne ) , Y. Ellinger Laboratoire d’Etude Theorique des Milieux Extremes, Ecole Normale Superieure, 24 rue Lhomond, F-75231 Paris Cedex 05, France ´ ˆ ´ Received 7 December 1998; in final form 6 April 1999
Abstract The structure of van der Waals dimers ŽN2 . 2 is studied using ab initio and density functional calculations. The potential energy surfaces corrected a priori for basis set superposition errors are necessary for determining both the geometries and the vibrational frequencies. With both ab initio MP2, MP4 and DFT PW91–PW91 levels of theory, the T-shaped and canted conformations appear to be the most stable, within 1–5 cmy1 of each other. The DFT PW91–PW91 dissociation energy is 67 cmy1 and an upper limit to the barrier to internal motion is 30 cmy1, both in excellent agreement with the values deduced from IR measurements. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction A large part of the mass of the universe is concentrated in dust particles. Due to the low temperature Ž10–50 K. of a number of celestial objects Žcold molecular clouds, comets, rings of giant planets, etc.. these particles behave as cold surfaces on which light molecules condense, forming what is known as ice mantles. This generic term, besides water, covers all types of organic or non-organic materials such as CH 4 , NH 3 , CH 3 OH, CO, CO 2 , or N2 whose accretion at surfaces is at the origin of the growth of interstellar grains. Before dealing with the structure of molecular solids and the possible chemistry that takes place at their surface under the extreme conditions prevailing in space, we found it necessary to consider the primary process of accretion, namely the formation of isolated dimers. )
Corresponding author. Present address: Department of Physical Chemistry, University of Geneva, 30 quai Ansermet, CH-1211 Geneva 4, Switzerland. E-mail: [email protected]
The first dimer, the study of which is reported in this Letter, is ŽN2 . 2 ; it may play an important role in the atmospheric chemistry of Neptune and Titan. Experimental evidence is limited but suggests that the most stable dimer has either a T-shaped or X-shaped geometry, with the centers of mass of the ˚ monomers separated by approximately 3.7–4.2 A w1,2x. From a theoretical point of view, the ŽN2 . 2 dimers have been studied with ab initio w3–8x and model potential w9–12x calculations. Most of these studies have predicted up to four local minima ŽTshaped, X-shaped, parallel and canted. but the energetic ordering varies with the method. The Møller–Plesset level of theory is usually accepted as a good starting point to account for the correlation effects required to describe weak interactions in van der Waals systems w13–15x. We will use the results obtained with this methodology for comparison with DFT calculations on this type of compounds. Although DFT gives remarkable results for solid structures as well as for covalent and ionic systems, it is not completely clear yet if it can
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 4 3 1 - 5
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
72
systematically deal with van der Waals supermolecules w16,17x even if clear-cut results have been obtained for rare gas dimers w18–20x, complexes between molecular nitrogen and benzene w21x and Cl 2 –C 2 H 4 w20x. Beyond the test case of the ŽN2 . 2 dimer, the present work is intended to probe the quality of DFT calculations in view of future studies on N2 , CO and CO 2 ices.
energy of the dimer in the molecular basis BD at geometry r D . The second term is a summation running over the energies of the monomers, each being calculated in the basis BD at their geometries in the dimer. The last term is a summation running over the energies of the monomers, calculated in the basis BM of each monomer at their geometries in the dimer. The BSSE-corrected PES is used to determine the geometry as well as the vibrational frequencies. The binding energy is given by:
2. Computational background D E s Ecorr y The ab initio calculations were performed at the MP2 and MP4 levels of theory as implemented in the GAUSSIAN94 package w22x. The standard 6-31 q G ) basis set was used. It has already proven reliable for these types of problems as well as for the description of molecular systems with a diffuse wavefunction such as negative ions w23x. Density functional theory ŽDFT. calculations were performed applying the Kohn–Sham ŽKS. formalism w24x. Two parametrizations of the exchange-correlation functionals which include gradient corrections within the generalized gradient approximation ŽGGA. were used: B3LYP w25,26x and PW91–PW91 w27x. DFT calculations were performed with GAUSSIAN98 w28x. In the DFT calculations, the atomic orbitals were constructed from the following contraction patterns: Ž7111r411r1) . w29x. This basis set is the closest to the basis used in ab initio calculations as shown by the comparison made with the results previously obtained on weakly interacting species w21x. The grids considered for the computations comprised 128 radial shells. At the end of the SCF procedure, each radial shell had 302 angular points. One of the main characteristics of the present study is that the considered potential energy surface ŽPES. is corrected for the basis set superposition error ŽBSSE. according to the counterpoise method w30x. The corrected energy is obtained as: Ecorr s ED Ž BD ,r D . y
Ý Ž EM Ž BD ,r D . Ms1,2
yEM Ž BM ,r D . . ,
Ž 1.
where the indices M and D refer to the monomer and dimer, respectively. The first term stands for the
Ý
EM Ž BM ,r M . .
Ž 2.
Ms1,2
It is now well established that this correction is essential for the determination of molecular structures bound by van der Waals interactions w31x. But unlike many studies where the BSSE correction is applied to the equilibrium structure, i.e. after the minimization has been completed, it is applied here to the entire optimization process. Therefore, we obtain a different surface, the minimum of which is the correct equilibrium geometry. We will show that, in these compounds where cohesion is dependent on weak intermolecular forces, neglecting BSSE or applying an inadequate correction can lead to dramatically erroneous results. As the BSSE-corrected energy is a sum of energies, it is possible to define a corrected Hessian H corr as the matrix of the second derivatives of the BSSEcorrected energy with respect to the nuclear coordinates, namely w32x: H corr s H D Ž BD ,r D . y
Ý ŽH M Ž BD ,r D . Ms1,2
yH M Ž BM ,r D . . ,
Ž 3.
where the terminology is the same as for the corrected energy and H indicates a Hessian matrix. Each second derivative matrix can be computed either analytically or numerically. Frequencies and normal modes corrected for BSSE are then obtained by diagonalization of the force constant matrix Fcorr which is derived from H corr . In such a way, the zero-point energy ŽZPE. obtained from these frequencies is also corrected for BSSE.
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
73
3. Results The geometry of the dimer in any conformation investigated can be defined by four parameters, namely the distance between the centers of mass R, two angles u 1 and u 2 , and one dihedral angle w ŽFig. 1.. Six types of conformations of ŽN2 . 2 have been studied: linear, parallel, T-shaped, X-shaped, canted and crystal ŽFig. 2., constrained to D`h , C2 Õ , D2 h , D2 d and C1 symmetries, respectively. 3.1. Symmetry constrained optimization The first step in our study consists in an optimization along the R coordinate. In these calculations the internuclear distance in the N2 molecule is the same for the dimer as for the monomer. The position of the minimum R 0 and the binding energy Emin are obtained by fitting n-order Ž n s 2–10. polynomials to 3–11 points. The BSSE correction appears to be essential for a reliable evaluation of both the equilibrium geometry and ZPE ŽFig. 4a.. More precisely, Table 1 shows that R min is systematically shifted toward longer distances, with a ˚ according to the correction varying from 0.2 to 0.4 A type of dimer. The smaller corrections appear in the DFT treatment. It can also be seen that the binding energy is systematically larger with the PW91–PW91 functional than in ab initio MP2 and MP4 calculations. In all cases it is greatly reduced when the BSSE correction is applied Žabout 80% for all conformations.. The BSSE-corrected potential energy curves obtained using MP2, MP4, DFT B3LYP and DFT PW91–PW91 are presented in Fig. 3 for the T-shaped conformation which will be shown below
Fig. 1. General representation for a dimer. M 1 and M 2 are the middle points of monomers 1 and 2. R is the distance between these two points. u 1 and u 2 are the angles between ŽM 1 M 2 . and the molecular axis. w is the dihedral angle between the two monomers.
Fig. 2. Different types of ŽN2 . 2 structures studied, as defined by the triplet Ž u 1 , u 2 , w ..
to be one of the two stable structures of equal energies. At the MP2 level, all dimers are found to be local minima, the T-shaped and canted conformations being the most stable ones with equivalent energies. Parallel and X-shaped conformations are found to be very close to each other, and the linear conformation Table 1 Equilibrium distances and binding energies along the R coordi˚. nate Ž D E in cmy1 and R min in A BSSE uncorrected
BSSE corrected
R min
DE
R min
DŽEqBSSE.
4.09 4.05 3.55 3.70 4.87
273 225 248 203 154
4.33 4.29 3.97 3.99 5.19
69 70 39 36 15
4.11 4.11 3.64 3.73 4.84
266 217 221 191 181
4.40 4.35 4.11 4.10 5.17
58 59 28 28 16
3.99 3.90 3.50 3.57 4.80
305 296 264 234 195
4.25 4.17 3.81 3.86 4.94
132 138 103 88 16
)
MP2r6-31qG T-shaped Canted X-shaped Parallel Linear MP4r6-31qG ) T-shaped Canted X-shaped Parallel Linear PW91–PW91r Ž7111r411r1) . T-shaped Canted X-shaped Parallel Linear
74
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
the dimer and will be discussed later when comparing larger N2 clusters to the crystal structure. 3.2. Minimum energy structures and Õibrational frequencies
Fig. 3. BSSE-corrected potential energy curves along the intermolecular coordinate R Žsee Fig. 1. for the T-shaped structure: MP2 Ž^., MP4 Žq., B3LYP Že. and PW91–PW91 Ž=..
is the least stable one. Binding energies at the MP4 level slightly differ from MP2, but do not interchange the relative stabilities of the dimers. For DFT calculations, all dimers were found to be non-bonded at the B3LYP level, implying that this functional is unable to reproduce dispersion forces for this type of molecular dimer. On the other hand, the PW91–PW91 results compare much better with the MP2 ones. The positions of the minima are close, ˚ toalthough systematically shifted by 0.10–0.15 A ˚ wards shorter distances Žexcept the 0.25 A shift for the linear dimer., the binding energy being a little larger. Very similar results have previously been obtained for the N2-benzene van der Waals complex w21x. The crystal conformation, as it appears in the a phase ŽFig. 2. is the third most stable arrangement after the T-shaped and canted conformations Ž D E s 56 cmy1 at the MP2 q BSSE level.. The distance between the two centers of mass optimized at the ˚ The MP2 q BSSE level of theory is R min s 4.2 A. experimental distance in the a phase of the crystal is ˚ Žcorresponding to a lattice parameter a s 5.7 4.0 A ˚A.. The a phase crystal conformation Ž u 1 s 90, u 2 s 35.3, w s 54.7. is not far from either the T-shaped Ž u 1 s 0, u 2 s 90, w s 0. or canted Ž u 1 s 45, u 2 s 45, w s 0. conformation. It can then be inferred that molecules tend to adopt T-shaped andror canted conformations but twist due to crystal constraints. This situation, however, has no physical meaning for
In order to assess the stability of the different dimers, frequencies have been computed on the corresponding MP2, MP4 and PW91–PW91 BSSE-corrected PES. As seen above, the force constant matrix corrected for BSSE is derived from a sum of Hessian matrices. Each Hessian can be computed either analytically or numerically on the MP2 q BSSE and PW91–PW91q BSSE surfaces, but only numerically for the MP4 q BSSE one. For the numerical calculations, we used a two-point differentiation. In order to obtain reliable results, especially for low frequencies, we had to use small displacements Ž d x . together with high accuracy in the energy value Ž d Eaccuracy .. We found, from comparison between analytical and numerical computations at the MP2 ˚ and level, that the parameters d x s 0.005 A d Eaccuracy s 10y7 a.u. give the best match between the two methods. MP4 q BSSE frequencies were then computed with these parameters. Frequencies have been obtained for the linear, parallel, X-shaped, canted and T-shaped conformations. All the frequencies of the T-shaped and the canted conformations are real, which confirms that these structures are two real minima ŽTable 2.. By Table 2 Internal modes Žcmy1 . of the potential surface for the T-shaped and the canted structures. Geometries are equilibrium geometries for each surface and force constant matrices were computed on the PESs corrected a priori for BSSE. The MP2 and PW91–PW91 force constant matrices were computed analytically, and the MP4 force constant matrix were computed numerically Žsee text for details.
T-shaped Bending Stretching Bending Out of plane Canted Bending Stretching Bending Out of plane
n MP 2
n MP4
n PW91 – PW91
b2 a1 b2 b1
;5 ;5 26 26
;5 ;5 16 21
26 27 37 39
b2 a1 b2 b1
;5 ;5 10 40
;5 ;5 8 32
22 24 27 61
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77 Table 3 Comparison between theoretical binding energies for the T-shaped and the canted conformations Ženergies in cmy1 .
T-shaped DŽEqBSSE. DŽEqBSSE.qZPE corr Canted DŽEqBSSE. DŽEqBSSE.qZPE corr Experiment Ž De . w1x
MP2
MP4
PW91–PW91
69 38
58 34
132 67
70 40
59 34
138 71 66
contrast, all the angular frequencies for the parallel, X-shaped, linear conformations are imaginary, which implies that these conformations are not stable.
75
Table 2 shows the vibrational frequencies for the internal modes of the most stable conformations: T-shaped and canted. ZPEs can be evaluated from these frequencies. The ZPE is a non-negligible quantity when compared to the binding energy: it represents up to 50%. The final binding energies are given in Table 3. Particularly good agreement is obtained between the DFT PW91–PW91 value and that of 66 cmy1 deduced from IR measurements by Long et al. w1x. 3.3. A priori Õersus a posteriori BSSE correction of binding energy One of the main points in this study is the a priori inclusion of the BSSE correction in the minimization process. This approach is seldom used, which is entirely justified when large energy differences are implied. In the present study, the energy differences involved are too small for this technique to be used. Using the T-shaped dimer as an example ŽFig. 4a. shows the result of an a posteriori BSSE correction on a minimum corresponding to an uncorrected PES. At the MP2r6-31q G ) level reported here, the error in the dissociation energy is 13 cmy1 , which represents 20% of the binding energy. The equilibrium distance between the centers of mass is under˚ In other cases, such as the estimated by 0.24 A.
Table 4 Comparison between internal modes Žcmy1 . before and after BSSE correction of the potential surface for the T-shaped structure ŽD 2 h .. Geometries are equilibrium geometries for each surface. The MP2 and PW91–PW91 force constant matrices are computed analytically, and the MP4 force constant matrix are computed numerically Žsee text for details.
Fig. 4. Stretching of the T-shaped conformation at the MP2r6-31 qG ) level. The dashed line represents the uncorrected energy and the plain line the energy corrected for BSSE. Ža. Difference between a posteriori correction of the BSSE Ž2. on the minimum of the non-corrected surface Ž1. and a minimization process on the BSSE-corrected surface Ž=.. Žb. ZPE corrections on the BSSEuncorrected and BSSE-corrected PES.
Before BSSE correction of the PES Bending Stretching Bending Out of plane After BSSE correction of the PES Bending Stretching Bending Out of plane
n MP 2
n MP4
n PW91yPW91
b2 a1 b2 b1
23 44 50 63
34 44 50 64
33 34 58 64
b2 a1 b2 b1
;5 ;5 26 26
;5 ;5 16 21
26 27 37 39
76
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
X-shaped conformation, the error in R increases up ˚ to 0.4 A. Using non-coherent ZPE is even worse. Fig. 4b shows two types of vibrational corrections, one computed on the corrected potential surface and hereafter referred to as ZPE corr Žtop., the other obtained from the uncorrected surface Žbottom.. Significant differences appear between the frequencies depending on whether they are corrected for BSSE or not ŽTable 4.. After correction, the low frequencies corresponding to the motion of one N2 molecule with respect to the other are lower than before correction. It was found that the stretching of the N2 bond is not affected by the formation of the dimer Žas anticipated from the conservation of the N–N bond length when going from the isolated monomer to the dimerized species.. It can be inferred from the curvature of the potential surfaces on Fig. 4a, that the second derivatives, and thus, the force constants, decrease when the BSSE correction is taken into account. By definition, BSSE is always negative. Since its magnitude decreases monotonically with intermolecular separation, it is quite obvious that correcting for BSSE results in a decrease of the PES curvature. A shallower surface implies lower frequencies for intermolecular modes and consequently a lower ZPE correction. For instance, at the MP2r6-31q G ) level, ZPE corr amounts to 31 cmy1 in the T-shaped dimer whereas the value directly obtained from the uncorrected PES is equal to 88 cmy1 which is 57 cmy1 higher. The present study clearly shows that erroneous results are obtained in such cases due to a critical overestimation of the vibrational frequencies. Table 5 Binding energies for the T-shaped conformation: comparison between a priori and a posteriori correction of the BSSE Ženergies in cmy1 ; =: instable.
DE a posteriori DEqBSSE DEqBSSEqZPE a priori DŽEqBSSE. DŽEqBSSE.qZPE corr
MP2
MP4
PW91– PW91
273
266
305
56 =
41 =
103 9
69 38
59 34
132 67
It can be seen in Table 5 that the T-shaped dimer is no longer bound at the MP2 and MP4 levels in these conditions Žthe same situation is found for the canted conformation.. 4. Concluding remarks In this contribution, we report a study of the structure and stability of the ŽN2 . 2 dimer at canted, T-shaped, X-shaped, parallel and linear structures plus that taken from the crystal. We have compared ab initio MP2 and MP4 results to the DFT ones. The present work provides further confirmation that the B3LYP functional is unable to properly describe the bonding forces in weakly interacting systems. By contrast, the exchange-correlation functional PW91– PW91 gives results which are not only close to the MP values but, even more satisfactorily, are in excellent agreement with IR measurements. The PW91–PW91 exchange-correlation functional, although still being semi-local and not correctly describing the Ry6 behaviour of the interaction energy at large distances, provides a very good approximation in the range of medium distances characteristic of the equilibrium structures of weakly interacting molecules. Moreover, this rather good performance of the PW91–PW91 functional, compared to that of B3LYP, can be attributed to its exchange part w21x. We have shown that PESs and ZPEs should be determined with explicit consideration of BSSE in the optimization procedure and calculations of the vibrational frequencies, a posteriori BSSE corrections possibly lead to erroneous equilibrium geometries, ZPE and binding energies. Under these conditions, the T-shaped and canted dimers are the only stable structures at the MP2, MP4 and DFT PW91– PW91 levels of theory, contrary to previous studies where equilibrium geometries were determined using PESs which were not BSSE corrected. It should be mentioned at this point that the same two conformations were found separated by 5 cmy1 in a recent study at the CCSDŽT. level with a posteriori BSSE corrections w33x. With proper a priori correction, the binding energy of the ŽN2 . 2 dimer is 38 cmy1 at the MP2 and MP4 levels of theory and 67 cmy1 for the DFT PW91–PW91 functional, in remarkable agreement with experiment.
O. Couronne, Y. Ellingerr Chemical Physics Letters 306 (1999) 71–77
The energy difference between the two stable conformations Ž1–6 cmy1 all corrections considered. is beyond the accuracy which can reasonably be expected when internal motions are treated at the harmonic level Žno coupling between vibrations.. No significant barrier could be found between the Tshaped and canted conformations, in agreement with the CCSDŽT. study w33x, which suggests that both conformations belong to the same potential well of the internal motions. The next saddle points, Xshaped and parallel can be seen as transition points for out-of-plane and in-plane Žnon-concerted. motions. We have verified that there is a direct path of decreasing energy connecting the parallel and Xshaped structures to the minima. These structures are about 30 cmy1 Ž43 K. higher than the equilibrium minima at all levels of theory. This value is in agreement with the 15–30 cmy1 Ž22–43 K. deduced from IR experiments w1x and is consistent with the permanent flip of the molecules in the crystal. Acknowledgements The authors would like to thank T.A. Wesołowski and Prof. N.C. Handy for stimulating remarks and suggestions in the course of this study. The financial support of CNRS Programme National de Planetolo´ gie ŽPNP. is greatly acknowledged. Part of the calculations were performed at IDRIS and CNUSC computing facilities. References w1x C.A. Long, G. Henderson, G.E. Ewing, Chem. Phys. 2 Ž1973. 485. w2x F. Carnovale, J.B. Peel, R.G. Rothwell, J. Chem. Phys. 88 Ž1988. 642. w3x F.H. Ree, N.W. Winter, J. Chem. Phys. 73 Ž1980. 322. w4x J. Tennyson, A. van der Avoird, J. Chem. Phys. 77 Ž1982. 5664. w5x P.J. Hay, R.T. Pack, R.L. Martin, J. Chem. Phys. 81 Ž1984. 1360. w6x H.-J. Bohm, R. Ahlrichs, Mol. Phys. 55 Ž1985. 1159. ¨ w7x F. Uhlik, Z. Slanina, A. Hinchliffe, J. Mol. Struct. ŽTHEOCHEM. 282 Ž1993. 271. w8x S. Raugei, G. Cardini, V. Schettino, Mol. Phys. 95 Ž1998. 477. w9x B. Jonsson, G. Karistrom, ¨ ¨ J. Chem. Phys. 74 Ž1981. 2896. w10x M.S.H. Ling, M. Rigby, Mol. Phys. 51 Ž1984. 855. w11x G. Brocks, A. van der Avoird, Mol. Phys. 55 Ž1985. 11.
77
w12x K.A. Franken, C.E. Dykstra, J. Phys. Chem. 97 Ž1993. 11408. w13x G. Granucci, M. Persico, Chem. Phys. Lett. 205 Ž1993. 331. w14x P. Hobza, H.L. Selze, E.W. Schlag, Chem. Rev. 94 Ž1994. 1767. w15x G. Chalasinski, M. Szczesniak, Chem. Rev. 94 Ž1994. 1723. ´ ´ w16x J.M. Perez-Jorda, ´ ´ A.D. Becke, Chem. Phys. Lett. 233 Ž1995. 134. w17x E. Ruiz, D.R. Salahub, A. Vela, J. Phys. Chem. 100 Ž1996. 12265. w18x Y. Zhang, W. Pan, W. Yang, J. Chem. Phys. 107 Ž1997. 7921. w19x D.C. Patton, D.V. Porezag, M.R. Pederson, Phys. Rev. B 55 Ž1997. 7454. w20x C. Adamo, V. Barone, J. Chem. Phys. 108 Ž1998. 664. w21x T.A. Wesołowski, O. Parisel, Y. Ellinger, J. Weber, J. Phys. Chem. A 101 Ž1997. 7818. w22x M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheeseman, T. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. Al-Laham, 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, Revision C.3, Gaussian Inc., Pittsburgh, PA, 1995. w23x O. Parisel, Y. Ellinger, C. Giessner-Prettre, Chem. Phys. Lett. 250 Ž1996. 178. w24x W. Kohn, L.J. Sham, Phys. Rev. A 140 Ž1965. 1133. w25x A.D. Becke, J. Chem. Phys. 98 Ž1993. 5648. w26x C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 Ž1988. 785. w27x J.P. Perdew, Y. Wang, in: P. Ziesche, H. Eschrig ŽEds.., Electronic Structure of Solids’ 91, Academie Verlag, Berlin, 1991. w28x M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, GAUSSIAN 98, Revision A.1, Gaussian Inc., Pittsburgh, PA, 1998. w29x F. Sim, D.R. Salahub, S. Chin, M.J. Dupuis, J. Chem. Phys. 95 Ž1990. 4317. w30x S.F. Boys, F. Bernardi, Mol. Phys. 19 Ž1970. 553. w31x P. Hobza, R. Zahradnık, ´ Chem. Rev. 88 Ž1988. 871. w32x S. Simon, M.D.J.J. Dannenberg, J. Chem. Phys. 105 Ž1996. 11024. w33x A. Wada, H. Kanamori, S. Iwata, J. Chem. Phys. 109 Ž1998. 9434.