Gas-phase proton affinity of ozone: a computational test of the experimental mechanism

Gas-phase proton affinity of ozone: a computational test of the experimental mechanism

Journal of Molecular Structure (Theochem) 543 (2001) 115±122 www.elsevier.nl/locate/theochem Gas-phase proton af®nity of ozone: a computational test...

106KB Sizes 0 Downloads 21 Views

Journal of Molecular Structure (Theochem) 543 (2001) 115±122

www.elsevier.nl/locate/theochem

Gas-phase proton af®nity of ozone: a computational test of the experimental mechanism M. Ceotto, F.A. Gianturco* Dipartimento di Chimica, Citta Universitaria, 00185 Rome, Italy Received 10 November 2000; accepted 23 November 2000

Abstract Extended ab initio calculations are carried out for the isolated ozone molecule, in its ground electronic state, and for its protonated adduct. For both systems, the corresponding nuclear geometries of lowest total energies are obtained by multidimensional optimization and different levels of basis set quality and correlation correction are examined in the two cases. The ®nal results provide very good accord between gas-phase protonation energy data and calculated quantities. The extensive numerical experiments clearly underline the need for balanced correlation treatments in systems where both static and dynamic correlation effects are important. They further indicate that the existing kinetics experiments on the gas-phase proton af®nity of O3 are not likely to involve the electron transfer channel via a sampling of the conical intersection con®guration space that exists in this system. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Gas phase; Proton af®nity; Ozone molecule

1. Introduction Gas-phase protonation processes have been of interest for many years to experimental and theoretical research groups both because of their relevance in many organic reactions and for their apparent simplicity that should make them amenable to accurate theoretical evaluations. In particular, the (O3H) 1 has been considered to be a relevant and active species during the oxygenation of alkenes, in acid medium, by the action of O3 [1±3] and several experiments have also postulated its existence as a reactive intermediate in the acid-catalyzed ozonolysis of alkanes [4,5]. Since this protonated adduct could be considered to be the conjugated acid of one of the * Corresponding author. Tel.: 1390-6-499-13306; fax: 1390-6499-13305. E-mail address: [email protected] (F.A. Gianturco).

least basic simple molecules, there have been several attempts at isolating (and characterizing) such an elusive and important species. Cacace and Speranza, for instance, have detected the protonated ion via a mass spectrometric experiment [5,6], which employed the delicate technique of the Fourier transform ion cyclotron resonance and managed to come up with a proton af®nity (PA) of the ozone partner in the gas phase of 148 ^ 3 kcal mol 21 at a temperature of 298 K. They further found [7] the relative basicity of O3 to be such that the corresponding protonated complex, (O3H) 1, was more stable than COH 1 but less stable than COSH 1, a relative ordering of strength which differs from the reversed order found in the condensed phase [8]. The corresponding computational studies have involved several types of calculations and the PA values found with them ranged from 124 to 157 kcal mol 21 depending on their level of sophistication,

0166-1280/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0166-128 0(00)00850-2

116

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

which went from simple Hartree±Fock (HF) computations to DFT model evaluations and to con®guration interaction (CI) studies [9±11]. As no ®rm consensus was reached by these earlier studies on the actual PA value or on the most stable geometry of the protonated complex, we have recently completed [12] a multicon®guration-self-consistent-®eld (MCSCF) treatment which de®nitely found the most stable structure to be the planar con®guration with the proton in an open-chain form located on one of the wing oxygen atoms. Furthermore, we underlined the existence of speci®c regions of the physical potential energy surface (PES) with strong non-adiabatic couplings between different electronic states and where possible electron transfer (ET) effects could occur leading to the complex fragmentations into a molecular ion and neutral hydrogen. Thus, kinetics measurements of PA could detect the re-crossing channel, if occurring, and would therefore obtain for the latter datum a lower energy value [12]. We have further carried out a more detailed analysis of such non-adiabatic features by studying the presence and the electronic nature of distinct regions of a reduced parametric space where conical intersections (CoIn) were located for at least two of the protonated adduct forms [13]. The experimental analysis [5], on the other hand, in general makes no mention of such electronic non-adiabatic features and considers the protonation process as the direct result of the reaction: O3 …gas† 1 H1 …gas† , O3 H1 …gas†

…1†

It therefore follows from the above that the question of obtaining a very accurate evaluation of the PA for the title system under computational conditions, which would allow a direct comparison with experiments has not been really answered in previous work and therefore it is still of interest to revisit the ionic complex of ozone with the aim of analysing as accurately as possible the computational feasibility of producing, for this four-atom system, the best available gas phase PA according to Eq. (1). The result could then be compared, in accuracy and in absolute value, with the existing experimental data mentioned before [5±7]. If in them the basicity kinetics, with respect to a reference species B, were to proceed via the break-up of the initial protonated adduct by also populating the ET channel through the physical

sampling of the CoIn con®guration space, then one could write the latter process via the following sequence, where the mechanisms invoked appear over the relevant arrows [12] and where the acronym CT refers now to the ®nal charge-transfer process to the reference base compound B used in the experiments: ET

¼HŠ 1 ‰BŠ ‰O3 HŠ1 1 ‰BŠ ! ‰O1 3 CoIn

CT

¼H±BŠ O O3 ¼‰H 2 BŠ1 ! O3 1 ‰HBŠ1 ! ‰O1 3 …2† The calculations with the direct mechanism of Eq. (1) should be able to show, by disagreeing with measurements and by producing values larger than theirs, that Eq. (2) indeed describes the correct physical process. The theoretical treatment is brie¯y outlined in Section 2, while Section 3 reports and analyses our present calculations. 2. The computational schemes The ozone molecule is known to be a computationally dif®cult system due to its diradical character, shown by the spatial shape of its 1a1 molecular orbital (MO), nearly exclusively located on the two wing atoms, and by the further existence of an equally important, alternative electronic con®guration in which the 1b1 MO is occupied instead. Hence, the effects of both static and dynamic electronic correlation need to be considered in a carefully balanced way for the initial neutral molecule [14]. Furthermore, when one examines the protonated adduct, we already found [12], in agreement with earlier studies [10], that the protonation process preferentially stabilizes the open-chain form with respect to the ring form (with an energy difference of about 23 kcal mol 21). Hence, the addition of the H 1 on the terminal, negatively charged oxygen atom tends to localize the charge to that end of the molecule, thereby partially destroying the O3 diradical character. The whole electronic structure of the complex is now frozen into one of the two, nearly degenerate dominant structures of the molecular partner [10±13]. Therefore, in order to be able to describe as best as

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

117

Table 1 The O3H 1 PA using the cc-pVTZ basis set expansion Data

MCSCF (close)

DFT-B3LYP

MRCI

CCSD(T)

ETOT(O3) (a.u.) a

2224.56969286

2225.49907126

2224.87019694

2224.83411334

2224.13389580

2224.11723995

2224.12749029

2224.13497704

2273.47

2867.11

2446.05

2438.71

0

±

±

O3 (a.u.) EHF

O3 b (kcal mol 21) Ecorr

CP correction (a.u.)

10.0006693

(kcal mol 21) ETOT(O3H 1) (a.u.)

10.42 2224.84691415

1

O3 H (a.u.) EHF 1

O3 H (kcal mol 21) Ecorr:

CP corrected protonation (a.u.) Energy (PE) (kcal mol 21) c 5 2 RT

(kcal mol 21)

[Evib(O3H 1)2Evib(O3)]

0 2225.74500020

± 2225.11760722

± 2225.08183726

2224.60810151

2224.61408889

2224.61094191

2224.60994220

2149.86

2709.66

2317.94

2296.12

20.27655199

20.24592894

20.24741028

20.24772392

2173.54

2154.32

2155.25

2155.45

11.48

11.48

11.48

11.48

17.58 d

17.58 d

17.58 d

17.58 d

21

(kcal mol ) O3 H1 O3 ‰Ecorr: 2 Ecorr: Š (kcal mol 21)

123.61

157.45

128.11

142.59

2DH (kcal mol 21) ˆ PA e

164.48

145.26

146.19

146.39

TDS (kcal mol 21)

17.76 f

17.76 f

17.76 f

17.76 f

156.72

137.50

138.43

138.63

21

2DG (kcal mol ) ˆ GB a b c d e f g

g

Ozone optimized geometry. tot tot Ecorr: ˆ Ecorrelated 2 EHF for each, different correlated calculation. …5=2†RT ˆ D…PV† 1 …3=2†RT: From the full computation of the vibrational modes using the DFT calculations. PAexp: ˆ …148 ^ 3† kcal mol21 : Estimated with the formula DS ˆ R ln ‰…s…O 3 † 1 s…H 1 †=s…O 3 H1 ††Š: GBexp : ˆ …140 ^ 3† kcal mol21 :

possible the marked correlation energy changes occurring in the adduct with respect to the initial neutral molecule we were forced to employ a variety of con®guration interaction approaches and to test their capabilities for this speci®c problem. Speci®cally, we have also employed the density functional approach (DFT) using the exchangecorrelation prescription of the B3LYP-model [15] with a basis set expansion which was increased from the cc-pVTZ to the cc-pVQZ quality. In our computations at the MCSCF level, the orbital space was divided into an internal space of 10 orbitals and this internal space further contained a region of core orbitals, which were doubly occupied, kept frozen and uncorrected during the CI process. It could also contain a region of closed-shell orbitals, which were doubly occupied but included in the correlation process. Finally, a space of active orbitals, was

selected, which were employed to de®ne the active space for the CI procedure. In all our present calculations, the active space was extended to all 12 orbitals employed within the CASSCF procedure, using the second-order direct multicon®gurational method of the molpro suite of programs [16±18]. Our results will be labelled here, as before [13], the MCSCF (close) calculations. We additionally carried out multireference-CI (MRCI) calculations to evaluate explicitly the correlation contributions used to improve on the MCSCF results [18]. We additionally carried out more extensive calculations using the coupled-cluster (CC) approach that included singles, doubles and estimated triples [18], labelled here as the CCSD(T) results. As a general indication, the total number of uncontracted con®gurations employed by the MRCI results was about 181 M for the protonated adduct and

118

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

Table 2 The O3H 1 PA using the cc2pVQZ basis expansion Data

MCSCF (close)

DFT2B3LYP

ETOT(O3) (a.u.) a

2224.58608588

2225.519577290

2224.14962268

2224.13276983

2273.88

2870.23

O3 (a.u.) EHF

O3 b (kcal mol 21) Ecorr:

CP correction (a.u.)

10.0002550

(kcal mol 21) ETOT(O3H 1) (a.u.)

10.16 2224.86357810

1

O3 H (a.u.) EHF 1

O3 H (kcal mol 21) Ecorr:

CP corrected protonation (a.u.) Energy (PE) (kcal mol 21) c 5 2 RT

(kcal mol 21)

0 0 2225.764064969

2224.62510017

2224.63120795

2149.65

2710.88

20.27723722 2173.97

20.244487679 2153.42

11.48

11.48

17.56 d

17.56 d

O3 H O3 ‰Ecorr: 2 Ecorr: Š (kcal mol 21) 2DH (kcal mol 21) ˆ PA e

124.23 164.93

159.35 144.38

TDS (kcal mol 21)

17.76 f

17.76 f

157.17

136.62

[Evib(O3H 1)2Evib(O3)]1 (kcal mol 21) 1

21

2DG (kcal mol ) ˆ GB a b c d e f g

g

Ozone optimized geometry. tot tot Ecorr: ˆ Ecorrelated 2 EHF for each, different correlated calculation. …5=2†RT ˆ D…PV† 1 …3=2†RT: From full calculations of the vibrational modes using the DFT calculations. PAexp: ˆ …148 ^ 3† kcal mol21 : Estimated with the formula DS ˆ R ln‰…s…O 3 † 1 s…H 1 †=s…O 3 H1 ††Š: GBexp: ˆ …140 ^ 3† kcal mol21 :

about 19 M for the neutral system. The reason for such an extensive series of comparative calculations was given by: (i) the desire to clarify via the present study the crucial role of dynamical and static correlation contributions for both the species involved in the gas-phase PA of ozone; and (ii) provide evidence that one can realistically accept Eq. (1) as the direct protonation mechanism for the present system and that there is no experimental evidence for the participation in it of the CoIn and ET channels described in Eq. (2).

3. Analysis of the results In order to obtain data, which could be meaningfully compared with the available experiments, we need to include several further corrections into the evaluation of the theoretical PA. The results reported

in Tables 1 and 2 include all such corrections which we are therefore explaining below: (i) all values of total electronic energies, for both the neutral molecule and the protonated complex, were computed using the fully optimized geometries given by each level of computation considered; (ii) to better underline the correlation contributions afforded by each level of computation, the results at the HF level are also reported in the tables; (iii) the basis-set-superposition-error (BSSE) correction is included via the Counterpoise method (CP) as already discussed in detail in our earlier work [13]. In most of the cases examined, and as expected from the size of the basis sets employed, the CP corrections turned out to be rather negligible; (iv) the vibrational modes were evaluated within

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

the DFT procedure in order to estimate the changes in vibrational energy content in going from the neutral molecule to the protonated adduct. The conventional de®nition of the PA in the gas phase gives it as the negative of the enthalpy change during the direct protonation process [11] of Eq. (1) PA ˆ 2DH ˆ 2DV 2 D…PV† ˆ 2DE 2 ‰Eint …O3 H†1 2 Eint …O3 †Š 2

5 RT 2

119

methods Ð the relative variations within the bond distances and bond angles remain within 1% of each other and mostly well below it. Thus, we can quantitatively see from our present theoretical data that the diradical character of ozone is markedly modi®ed upon protonation into the openchain form and therefore dynamical correlation effects are important for it. 2. The above point is given a more quantitative quali®cation when we compare the correlation

…3†

where the last term on the r.h.s. of Eq. (3) is the sum of the D(PV) value and the energy difference in the proton translational energy contributions to internal energy, before and after the attachment in the gas phase. It mainly corresponds to the value of the proton translational energy content before attachment and is reported as such in Tables 1 and 2. We further estimated the entropy variation during the process by de®ning it as follows 2TDS ˆ 2T‰SO3 H1 2 SO3 2 SH1 Š

…4†

which we evaluated explicitly by considering all the entropy contributions for each of the optimized structures of the isolated molecule and for the protonated complex (see footnote to both Tables 1 and 2). We could therefore obtain, from the above estimates, the full value of the DG of the reaction of protonation, a quantity which de®nes the gas basicity (GB) of the molecular species under consideration DG ˆ GB ˆ PA 2 TDS

…5†

If we now analyse the results reported in Table 1, we can make the following observations: 1. All the four correlated computations also provide the optimized, Born±Oppenheimer (BO) geometries for the most stable form of the protonated complex Ð the open-chain proton attachment to a ring oxygen. In order to compare the details of such optimized geometries, we present them in Fig. 1, where the Ê ) and of the bond actual values of the distances (in A angles (in degrees) are also given. The data in the ®gure clearly show the remarkable agreement between the structural data provided by the four

Fig. 1. Computed geometries for the most stable nuclear con®gurations obtained for the open-chain form of the protonated adduct of ozone, using the computational procedures discussed in the present work. The acronyms are described in the main text. Distances are in Ê and angles in degrees. A

120

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

energy corrections which are afforded by the various levels of our calculations. They are shown in the third and seventh lines of Tables 1 and 2. The DFT results provide dramatically larger values for such correction, indicating the wellknown feature of the method that consistently overestimates correlation corrections. On the other hand, the correlation corrections for the neutral molecule turn out to be, with all methods, larger than those obtained for the protonated complex, thus con®rming that the disappearance of the dicationic character in going to the open-chain form of the adduct somewhat reduces correlation effects by localizing paired electrons onto either wing atoms. 3. The procedure which yields the least variation in the correlation correction is the MCSCF treatment. Its results also change very little when one goes from the cc-pVTZ basis set to the more extended cc-pVQZ basis set. As a consequence of this increased correlation rigidity upon protonation, we see that the energy contributions coming from correlation energy changes in going from the free species to the protonated complex turn out to be the smallest in the case of the MCSCF results and the largest for the DFT calculations (see the 11th lines in Tables 1 and 2). 4. The data reported in the Tables 1 and 2 also indicate that the MRCI correlation-corrected procedure and the more extended CCSD(T) calculations yield the lowest total energy for both the free ozone molecule and for the protonated adduct, although they start with HF energies for ozone which are very close to those yielded by the MCSCF (close) procedure. However, when we look at the total electronic energies which are obtained for the complex, we see that the MRCI and CCDS(T) calculations have already improved on the MCSCF (close) results at the HF level, becoming markedly better than them when correlation corrections are added. 5. It is also interesting to note that the BSSE correction included via the CP method [19], the CP energy corrections shown in Tables 1 and 2, is really negligible in the case of the free ozone molecule where the large basis sets employed saturate the system, while it is not negligible in the case of the ionic complex. The fact that the MCSCF

calculations require the largest CP correction is also indicative of its difference in quality when describing the correlation correction of the complex. 6. If we now examine the ®nal computed PA in Tables 1 and 2, we see that it is very strongly affected by the changes in correlation energy contributions that each computational procedure generates. Thus, the MCSCF procedure, which yields a reduced description of such changes, is not accurate enough to provide PA values within the experimental error, not even when a larger basis set is employed. On the other hand, both MRCI and CCSD(T) calculations show remarkably good agreement with experiments since they provide larger contributions from correlation energy changes and less marked CP corrections. It is also interesting to see how the DFT calculations, by compensation of effects, also yield a PA value in very close agreement to experiments. In all instances, however, the computation of the PA via Eq. (1) appears to be the correct description of the physical process. To further verify how effective the methods could be in evaluating the asymptotic energy gap between the two main dissociation channels from the ground electronic state of (O3H) 1, i.e. the O3 1 H 1 and the O1 3 1 H; both of which we have already discussed in our earlier work [12,13], we carried out additional calculations which can be compared with the experimental value for the energy separation. This quantity is given by the difference in ionization potentials between the molecular partners [12], 1.07 eV. The calculations reported in Table 3 show that, once again, the DFT method provides very good agreement with the experimental value, this being true for both the basis set choices which we have considered. On the other hand, the MCSCF calculations appear to worsen when the basis set is enlarged since the description of the ionic ozone molecule requires strong correlation changes which are not given well by the MC procedure. When we examine the outcomes from the MRCI and the CCSD(T) calculations, we see clearly their improved treatment of correlation effects, especially for the case of the larger basis sets, hence a better agreement with the experimental data.

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

121

Table 3 Calculation of asymptotic energy gap for the protonation reaction Data

MCSCF (close)-vtz

DFT-B3LYP-vtz

MRCI-vtz

ETOT(O3) (a.u.) a EO3 1 EH1 (a.u.) EO13 (a.u.) a EO13 1 EH (a.u.) …EO3 1 EH1 †2 …EO13 1 EH † (eV) b …EO3 1 EH1 †2 …EO13 1 EH † (kcal mol 21) c

2224.56969286 2224.56969286 2224.11225066 2224.61225066 11.15

2225.499071261 2225.499071261 2225.038980945 2225.538980945 11.09

2224.87019694 2224.87019694 2224.40540690 2224.90540690 10.96

Data ETOT(O3) (a.u.) a EO3 1 EH1 (a.u.) EO13 (a.u.) a EO13 1 EH (a.u.) …EO3 1 EH1 †2 …EO13 1 EH † (eV) b …EO3 1 EH1 †2 …EO13 1 EH † (kcal mol 21) c

MCSCF (close)-vqz 2224.58608588 2224.58608588 2224.17093616 2224.67093616 12.31

a b c

126.70

153.74

125.04 DFT-B3LYP-vqz 2225.519577290 2225.519577290 2225.057729042 2225.557729042 11.04 123.94

122.09 CCSD(T)-vqz 2224.864227034410 2224.864227034410 2224.400327227055 2224.900327227055 11.04 123.88

Optimized geometry. Eexp: ˆ 11:07 eV: Eexp: ˆ 124:639 kcal mol21 :

4. Conclusions In the present work we have analysed in some detail the calculation of the gas-phase PA for a special molecule like ozone. The experimental quantity is known but an independent theoretical con®rmation of its values is still lacking. Thus, we have carried out here a sort of computational experiment where different procedures are employed to treat the various effects which contribute to PA, especially the very important correlation changes induced in the ozone molecule via the charge localization created by the direct proton attachment. All the calculations reported here include a rather large portion of the correlation corrections, both in the free fragments and in the complex. Hence, the present study is able to show that the dicationic character of ozone undergoes marked modi®cations in the protonated adduct and therefore it requires a very delicate balance between different correlation effects in the two species involved in order to produce computed PA values which agree with the experimental ®ndings. Hence, we ®nd that fuller correlation procedures as those used in the MRCI and the

CCSD(T) calculations correctly include balanced correlation corrections in free ozone and in the protonated complex, thereby yielding remarkably good agreement with experiments. Such control of the quality of correlation changes is also shown to be needed when ionic ozone is produced in its ground electronic state in order to obtain the correct energy gap between the two ET and CT channels considered here [12]. The treatment of the protonation reaction as a direct attachment process as in Eq. (1), and therefore the exclusion of possible multistep contributions from the ET channel …O1 3 1 H†; as reported in the sequence of Eq. (3), when observing the break-up of the complex in a kinetic mass spectrometric experiment (also suggested to exist in other systems [20]) is clearly indicated by the present calculations. We have thus been able to reproduce the measured data in a reasoned way and to account in detail for the various physical effects that are likely to contribute to the ®nal value of the gas-phase PA of the present molecule. Our calculations further show that the kinetic process occurring during the experimental measurements can be taken to be a one-step transfer process to the reference base B and therefore we can exclude that

122

M. Ceotto, F.A. Gianturco / Journal of Molecular Structure (Theochem) 543 (2001) 115±122

the CoIn contributions to the dynamics play any signi®cant role in producing the ®nal PA value. Acknowledgements The ®nancial support of the Italian National Research Council (CNR), the Italian Ministry for University and Research (MURST), the Research Committee for the University of Rome and the Max±Planck Research Society is gratefully acknowledged. We are also grateful to professor Speranza and professor Cacace for several illuminating discussions on the experimental ®ndings. References [1] G.A. Olah, N. Yaneda, D.G. Parker, J. Am. Chem. Soc. 98 (1976) 5261. [2] P.A. Nangia, S.W. Benson, J. Am. Chem. Soc. 102 (1980) 3105. [3] N. Yaneda, G.A. Olah, J. Am. Chem. Soc. 99 (1997) 3133.

[4] R. Cipollini, M.E. Crestoni, S.J. Fornarini, J. Am. Chem. Soc. 119 (1997) 9499. [5] F. Cacace, M. Speranza, Science 265 (1994) 208. [6] F. Cacace, M. Speranza, Inorg. Chem. 35 (1996) 6140. [7] F. Cacace, M. Speranza, Chem. Lett. (1998) 419. [8] L. Mariey, J. Lamotte, P. Hoggan, J.C. Lavalley, K. Boulanine, A. Tsyganenko, Chem. Lett. (1997) 835. [9] M. Kausch, P.J. R.-Schleyer, J. Comput. Chem. 1 (1980) 94. [10] C. Meredith, G.E. Quelch, H.F. Schaefer III, J. Am. Chem. Soc. 113 (1991) 1187. [11] M. Toscano, N. Russo, J. Rubio, J. Chem. Soc., Faraday Trans. 92 (1996) 2681. [12] M. Ceotto, F.A. Gianturco, J. Phys. Chem. A 103 (1999) 9984. [13] M. Ceotto, F.A. Gianturco, J. Chem. Phys. 112 (2000) 2658. [14] K.B. Mathisen, O. Grapen, P.N. Skancke, U. Wahlgren, Acta Chem. Scand. A37 (1983) 817. [15] R.G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, Oxford, 1989. [16] H.J. Werner, P.J. Knowles, J. Chem. Phys. 89 (1988) 5803. [17] H.J. Werner, P.J. Knowles, Chem. Phys. Lett. 145 (1988) 514. [18] H.J. Werner, Adv. Chem. Phys. 69 (1987) 399. [19] S.F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553. [20] R. Contreras, V.S. Safant, J. Andres, P. Perez, A. Aizman, O. Tapia, Theor. Chem. Acc. 99 (1998) 60.