A theoretical analysis of the acid–base equilibria of hydroxylamine in aqueous solution

Chemical Physics Letters 490 (2010) 159–164

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

A theoretical analysis of the acid–base equilibria of hydroxylamine in aqueous solution M.I. Fernández, M. Canle, M.V. García, J.A. Santaballa * Chemical Reactivity and Photoreactivity Group, Department of Physical Chemistry and Chem. Eng. I, University of A Coruña, Rúa Alejandro de la Sota 1, E-15008 A Coruña, Spain

a r t i c l e

i n f o

Article history: Received 4 November 2009 In final form 17 March 2010 Available online 19 March 2010

a b s t r a c t Ammonia oxide (+NH3O), a zwitterionic tautomer of hydroxylamine (NH2OH), has been proposed to explain the high reactivity of NH2OH with phosphate esters. The key parameter for kinetically significant reactions through +NH3O is the value of the tautomeric constant KT = [+NH3O]/[NH2OH], which it has so far proved elusive to measure experimentally. Thermodynamic cycles using CBS-QB3 and CPCM computed data have been used to estimate microscopic and macroscopic acid–base equilibria involving NH2OH, the latter in close agreement to the experimental ones, but those involving anions. Calculated KT is 2.6  102, which indicates that although molecular NH2OH is favoured, there should be enough + NH3O in aqueous solution to be a kinetically active species. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Knowledge of the pKa values of ionizable groups is of particular interest for elucidating reaction mechanisms, especially those involving proton transfers, and also for interpreting the way substrates bind to enzymes. Usually, macroscopic pKa values are known, while microscopic constants, especially tautomerization constants KT, are not. Their accurate determination is relevant for complete understanding of chemical behaviour and biological activity, but this is usually an extremely difficult task from the experimental point of view [1]; in such cases computational chemistry has revealed as an useful tool to estimate microscopic pKa values [2]. Although the chemical structure of NH2OH is very simple, its reactivity is rather complex. It is an unique, ambident nucleophile showing a-effect. Generally speaking this term refers to the reactivity enhancement observed for nucleophiles with lone pair electrons on the atom adjacent to the nucleophilic center. Alkylation takes place on the nitrogen, whereas usually it is both acylated and phosphorylated on oxygen, indicating reactivity with harder electrophiles is favoured through oxygen [3–5]. NH2OH is also of particular mechanistic interest because it is known to react with phosphate esters through oxygen, rather than through its more basic nitrogen center [3,4]. The mechanism of nucleophilic attack by NH2OH through oxygen has been discussed a number of times, but always inconclusively [4]. There is still an open debate on the role of ammonia oxide (+NH3O), the zwitterionic tautomer of NH2OH, in this mechanism [3,4,6]. Although

* Corresponding author. E-mail addresses: [email protected] (M.I. Fernández), [email protected] (M. Canle), [email protected] (M.V. García), [email protected] (J.A. Santaballa). 0009-2614/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2010.03.047

the existence of ammonia oxide has been discussed as early as 1896 [7], direct evidence of its existence is ‘a major challenge to experimentalists working in the field of elusive molecules’ [8]. Noticeably, +NH3O is not a typical oxyanion, because on the nitrogen it bears up to three hydrogens adequately located to facilitate a thermodynamically favourable proton transfer. Recently Kirby et al. reported the first firm evidence for the existence of +NH3O in the condensed phase, in the shape of a crystal structure containing equal amounts of the neutral and the conjugate acid 7 [6]. The same authors quite recently have estimated a value of KT = 0.22 (298 K) based on structure–activity relationships [9]. To the best of our knowledge, this tautomeric equilibrium has so far proved impossible to be experimentally measured, so we decided to carry out theoretical calculations to estimate its equilibrium constant. Although NH2OH and some related species have been submitted to calculation [10–16], few works deal with the estimation of their tautomeric constants, the neutral form being by far the predominant species [6,17]. In this work in addition to the calculation of the tautomeric equilibrium between NH2OH and +NH3O a complete set of acid– base microscopic equilibrium constants, and the corresponding macroscopic ones, involving NH2OH are calculated by using the compound method CBS-QB3, to compute energies in the gas phase, and single point PCM calculations to include aqueous solvation. 2. Computational methods All calculations were carried out using the GAUSSIAN03 suite of programs [18]. Absolute pKa determinations are dependent on accuracy, therefore the choice of theory level is critical. Compound methods like Gx family or complete basis set (CBS) are usually the

160

M.I. Fernández et al. / Chemical Physics Letters 490 (2010) 159–164

methods of choice to obtain vacuum energies [19,20]. Here, CBSQB3 was selected as it has been the reference method in recent work on pKa calculations [21]. This method carries out an initial geometry optimization and frequency calculation using the B3LYP/CBSB7 level, followed by single-point calculations at the CCSD(T)/6-31+G(d), MP4SDQ/CBSB4, and MP2/CBSB3 levels. The final single-point calculation implies a complete basis set extrapolation. This method includes the zero-point vibrational energy (ZPVE) correction from the B3LYP/CBSB7 frequency calculation scaled by 0.9899. The geometries have been fully optimized, and the nature of the stationary points was checked by looking at the imaginary frequencies, zero for minima and one for first-order saddle points. Solvation free energies were calculated with the conductor-like polarizable continuum model (CPCM) [22], as implemented in GAUSSIAN03 [18]. Single-point energy calculations in solution, using the CBS-QB3 optimized geometries, were carried out at different levels of theory: CPCM/HF/6-31G(d), CPCM/HF/6-31+G(d,p), CPCM/HF/6-3111+G(d,p), CPCM/fullMP2/6-3111+G(d,p), CPCM/ B3LYP/6-31+G(d,p), and CPCM/B3LYP/6-3111+G(d,p). Average values were used in the calculation of pKas. As usual in CPM models, the molecular cavity involves the partial overlapping of spheres around each element, which are characterized by its radius, and are built as sum of small triangles called tesserae; their number and area were set at 240 and 0.3 (in Å2), respectively. UAO (default option in G03), UAKS and Pauling (actually Merz–Kollman) radii, as implemented in GAUSSIAN03 [18], were initially used.

3. Results and discussion 3.1. Energy and geometry In spite of the simplicity of the species involved, up to 6 atoms, some of them allow conformational isomers. Fig. 1 shows the main species and the conformational isomers considered in this work, and Table 1 contains the corresponding geometric parameters (see also the Supplementary material). NH2O (1) has a lower energy than HNOH (2) both in the gas phase and in solution (Table 2). Analysis in terms of the NBO theory indicates that the higher donor–acceptor delocalization takes place in 1, the strongest interaction taking place between one of the oxygen lone pairs, as donor, with N–H antibonding orbitals, as acceptors. On the other hand, taking NH2OH 4a as a reference, both the proton affinity (PA) and the gas phase basicity (GB) are lower, 6 and 10 kcal mol1 respectively, for NH2O than for  HNOH; which is consistent with the expected higher acidity of the O–H bond. The calculated N–O bond length in 1 is 1.430 Å (Table 1), typical of single N–O bonds, whereas the same bond in 2 is ca. 0.14 Å longer; the increase of the negative charge on the nitrogen atom in 2 leading to a weakening of this bond. Charges obtained from Mulliken and Natural Population Analysis (NPA) [23] agree with oxygen and nitrogen being the most negatively charged centers in 1 and 2, respectively (Table S1 of Supplementary material). At the CBS-QB3 level two rotational barriers have been found for 2, corresponding to syn and trans conformations, which lie 5.78 and 1.71 kcal mol1, respectively, above the minimum (see Supplementary material). NH2OH (4) has two conformational isomers, anti-NH2OH (4a) and syn-NH2OH (4b), both having the same total energy in solution: anti-NH2OH is 4.35 kcal mol1 more stable than syn-NH2OH in the gas phase, whereas the solvation energy of the latter is 4.34 kcal mol1 lower than that for the former (Table 2). Previous

Fig. 1. Structures and hydrogen atoms numbering of the hydroxylamine and its derivatives used in this work.

calculations lead to the same prediction in the gas phase, the anti-isomer being favoured [8,13]. Following the common convention, molecular NH2OH will be referred to as neutral species, although strictly speaking +NH3O is also neutral, with zero net charge as a result of two opposite charged ends (zwitterion). From the relevant geometric parameters of the optimized structures (Table 1), it follows that the calculated values for anti-NH2OH closely agree with the experimental values for NH2OH, obtained by microwave spectroscopy [24,25]. Calculated geometric parameters for 4a are similar to those obtained with other methods (Table 3); usually the main difference, not high, from the experimental values is the HNH angle. According to the NBO theory [23], the higher donor–acceptor delocalization takes place in syn-NH2OH, in both cases the strongest interaction being between one of the oxygen lone pairs, as donor, with N–H antibonding orbitals, as acceptors. In addition to conformational isomers, NH2OH has two struc+ H3NO (5), both zwitterions. tural isomers: HNOHþ 2 (3) and Two conformational isomers have been found for HNOHþ 2 (3): syn (3a) and anti (3b). As found in previous studies (Table 3), the

161

M.I. Fernández et al. / Chemical Physics Letters 490 (2010) 159–164 Table 1 Relevant geometric parameters of hydroxylamine and its conjugated acid and base derivatives. Level of calculation: CBS-QB3. dN–Oa

dN–Hb



NH2O 1

1.430



HNOH 2 HNOHþ 2 3a

1.574 1.679

1.045 1.045 1.042 1.038

HNOHþ 2 3b

1.718

1.031

NH2OH 4a

1.446 1.453d

1.019 1.019 1.016d

NH2OH 4b

1.430

0.967

+

1.354

1.019 1.019 1.039 1.039 1.039 1.023 1.022

0.978 0.980

107.7 107.7 114.0 114.0 114.0 105.3 101.8





a b c d

NH3O 5

dO–Hb

0.958 0.965 0.965 0.965 0.965 0.962 1.962d

107.9 107.9 97.2 95.1

103.8 103.8 103.3d

1.483

H2 NOHþ 2 6b

1.530

1.028 1.028

0.981 0.981

99.9 99.9

+

1.406

1.031 1.031 1.028

0.974

112.5 112.5 104.8

H3NOH 7

\HNOHc

\NOHc

\HNHc

\HOHc

98.6

93.7

H2 NOHþ 2 6a

106.8 59.7 60.0 126.4 126.0 124.9 57.9e 125.1 57.9d 57.9 57.9

98.1 109.9 109.8 96.9 96.9 102.1 101.3d

105.4 107.1d

107.9

107.6

98.8 106.5

104.5

24.7 105.8 90.0 139.5 69.8 176.4 176.8 69.8 62.6 62.6 180.0

118.5 118.8

109.9

112.8

106.9 106.9

104.7

106.6

107.2

110.3 108.2 108.2

In Å. Distances N–Hi (or O–Hi), in Å, ordered from 1 to n, as in Fig. 1. Angles \HiNOHj (\HiNHj and \HiOHj), in °, from the lowest value of i to the highest value of j according numbering depicted in Fig. 1. Exptl. value, Refs. [24,25].

Table 2 Computed thermodynamic data for hydroxylamine and its derivatives (reference compound for each group is underlined). Level of calculation: CBS-QB3 – gas phase, and CPCM – solvation. Species

D Ha (Hartrees)

DG a (Hartrees)

DGgasa (kcal mol1)

DGsolva (kcal mol1)

NH2O 1  NHOH 2  NHOHþ 2 3a

130.924309 130.907802 131.431645 131.438905

130.950288 130.934461 131.459337 131.466595

0 9.9 65.7 61.2

96.81 88.89 24.75 18.69

131.537437 131.530578 131.506573 131.802547

131.564074 131.557135 131.523508 131.829835

0 4.4 25.5 26.8

10.14 14.48 33.43 90.34

131.801531 131.845568

131.829810 131.872623

26.9 0

90.61 93.85

NHOHþ 2 3b NH2OH 4a NH2OH 4b NH3O 5 NH2 OHþ 2 6a NH2 OHþ 2 6b + NH3OH 7 

a

\HNOc

Species

T = 298.15 K.

N–O bond distance in 3 is abnormally large for a single N–O bond (>1.67 Å), 3b showing the longest N–O length (Table 1). In the gas phase the anti isomer (3b) is the most stable species by 4.55 kcal mol1, while the in aqueous solvation, the syn form (3a) becomes more stable by 1.50 kcal mol1 (Table 2). The energy difference between anti-NH2OH (4a) and both isomers (3a and 3b) is greater than 60 kcal mol1 (Table 2). Similar values have been obtained previously [8,13]. + H3NO (5) is the most symmetrical species (C3v), showing the shortest N–O bond distance (1.354 Å). In terms of the NBO theory, the higher donor–acceptor delocalization takes place between one of the oxygen lone pairs and the three N–H antibonding orbitals. This species lies 25 kcal mol1 above anti-NH2OH (4a) (Table 2), a value similar to others calculated previously [6,8,14,26]. When +H3NO protonates (7) the N–O bond distance increases to 1.406 Å, length still typical of a single N–O bond. In contrast to

+

H3NO, in its protonated form (7) the three N–H distances are not equal, the N–H bond, anti to the O–H bond, being shorter than the other two by 0.002 Å. Similar differences are obtained for other geometrical parameters involving those N–H bonds (Table 1). In terms of the NBO theory, the main donor–acceptor delocalization is observed between one of the oxygen lone pairs and the two equivalent N–H antibonding orbitals, which agrees with observing a difference between the three N–H distances. Structural isomers of 7 have also been calculated (6a and 6b). Both species, endergonic with respect to 7, have almost the same energy in the gas phase, and also when solvation is taken into account (Table 2). As expected the N–O bond length grows up, c.a. 1.48 Å (6a) and 1.53 Å (6b). Table 3 collects the computed N–O bond distances for some of the previous compounds, the differences arising from the level of calculation and basis set. All values are similar for the same compound, except for H2 NOHþ 2. Comparison between the energies of +NH3OH and H2 NOHþ 2 points to the preferred site of protonation being the N atom by c.a. 27 kcal mol1 (Table 2), an energy gap similar to previously calculated values [2,8,16,27]. 3.2. pKa calculation Theoretical pKa values have been calculated with the thermodynamic cycles shown in Fig. 2, for which the results must be consistent [28]. Although discussion could be open about the most rigorous way to theoretically estimate pKa values, from a practical point of view the challenge is to have reliable method(s) available to estimate them. One of the weak points of the current theoretical computation of pKa values is the calculation of accurate solvation energies for H+ or H3O+. The precision of the calculated pKa values partly relies on the accuracy of those values. The standard Gibbs free energy for

162

M.I. Fernández et al. / Chemical Physics Letters 490 (2010) 159–164

Table 3 Survey of computed nitrogen–oxygen bond distances (Å) for hydroxylamine and its derivatives from the literature. Ref.

HNOHþ 2 3a 

HNOHþ 2 3b 

[6] [8] [10] [11,12] [13] [14] [15] [16] [17]

1.632

1.66

1.69

NH2OH 4a 1.444 1.460 1.44 1.46

1.43

+

NH3O 5

1.437

1.367 1.419 1.352

1.44

1.69 1.377 1.36

1.4508

[26]

1.451 1.403

[31]

1.455 1.435 1.434 1.388 1.423 1.432 1.433

[32]

[33] [34] [35] [24,25] (exp.) This work

NH2OH 4b

H2 NOHþ 2 6a

H2 NOHþ 2 6b

1.418

1.479 1.5177

1.408 1.4106 1.4175 1.4083

1.441 1.461 1.482 1.486

1.367 1.397 1.406 1.407 1.374

1.362 1.376

1.395 1.434 1.444 1.444 1.403 1.4359

1.490 1.460

1.679

1.718

1.435 1.430

1.446

NH3OH+ 7

1.354

1.483

1.392 1.530

1.406

Fig. 2. Thermodynamic cycles used in the calculation of aqueous pKa values.

the formation of H+(g) was calculated by using the Sackur–Tetrode equation (6.28 kcal mol1), which reduces to 4.39 kcal mol1 when the change of reference state, from atmospheres to mol L1, is considered. The value for the solvation free energy of H+, still subject to debate, has been taken as –264 kcal mol1, and on this basis 110.8 kcal applies to H3O+ [29]. In Fig. 2b the solvation free energy for water was taken as 6.32 kcal mol1 [29]. Besides the mentioned problem in selecting proper solvation values, in the

past decade different pKa values have been calculated using the same solvation values [29]. Discrepancies are usually due to the use of different reference states [29]. Although 1 atm was chosen for the ideal gas standard state in this work, similar results would be obtained assuming a standard state in gas phase of 1 mol L1 [2]. Results obtained using the two cycles shown in Fig. 2 differ by only 0.10 units. In the remainder of this work reported values are

163

M.I. Fernández et al. / Chemical Physics Letters 490 (2010) 159–164 Table 4 Influence of solute cavity radii on the pKas calculation, and computed values for reference compounds. Level of calculation: CBS-QB3 – gas phase, and CPCM – solvation. Acid/base system

NHþ 4 /NH3 CH3 NHþ 3 /CH3NH2 NH3/NH 2 CH3NH2/CH3NH CH3OH/CH3O (CH3)3N+OH/(CH3)3N+O a

Exp.a

Radii UAO

UAKS

PAULING

1.55 5.47

5.55 6.25

10.24 13.14 35.14 40.12 25.43 4.90

9.25 10.66 35.00 35.00 15.50 4.65

Aqueous solution and T = 298 K.

the average of both cycles. Note that corrections should be taken into account when two or more species, say syn- and anti-NH2OH, are close enough in energy: the value used in the cycle should be an average value resulting from the Boltzmann’s distribution. On the other hand, no correction applies when the energy gap is high, or the species have the same total Gibbs free energy (4a and 4b, and 6a and 6b). In order to choose the best radius set, the pKas of NHþ 4 and CH3 NHþ 3 ions were calculated (Table 4), both the number of tesserae and their area having no influence on the results. From there it follows the best option is to use Pauling radii. Calculated values for our system at 298 K are shown using Pauling radii (see Fig. 3) and default radii (see Supplementary material), the latter being worse. As expected on chemical grounds, the relevant species participating in acid–base equilibria, both in aqueous solution and in the gas phase, are those in blue. The computed pKa (+NH3OH + H2O ¡ NH2OH + H3O+) is not far from the empirical value, whereas a noticeable difference is obtained when the deprotonation of hydroxylamine pKa (NH2OH + H2O ¡ NH2O + H3O+) is considered. pKas of reference compounds were estimated to check the reliability of this value (Table 4). Except in the case of NHþ 4, some improvement in the computation of solvation energies involving anions should be done to produce more accurate values. With this in mind, and to test the performance of the method when

Fig. 3. pKa and pKT values for hydroxylamine and its derivatives in aqueous solution at 298 K computed by using the composite CBS-QB3 method for gas phase energies, and CPCM to simulate aqueous solvation.

zwitterionic species are involved, pKa (+N(CH3)3OH + H2O ¡ (+N(CH3)3O + H3O+) was calculated, the obtained value being in good agreement with experiment (Table 4) [6]. pKas for the relevant acid–base equilibria of NH2OH in aqueous solution are closer to empirical ones than previously calculated [6,16] (see Table 5). The difference cannot be explained only by the use of the value 262.4 kcal mol1, instead of 264.0 kcal mol1, for the Gibbs free energy for proton salvation. The effect of such change is to make the calculated values ca. 1.2 pKa units more positive. Our computed values are also consistent with the corresponding pKas estimated from linear free energy relationships (LFER) (Table 5) [9]. Although computational results should be considered with caution, the value for the tautomeric equilibrium NH2OH ¡ +NH3O is

Table 5 Microscopic and macroscopic pKa values of hydroxylamine obtained by CBS-QB3 – gas phase, and CPCM – aqueous solvation.

Computational LFERb

pK1

pK3

5.21, 1.00 [17], 3.2 [17] 5.69 ± 0.45

20.73, 28.7 [17], 25.5 [17] 13.41 ± 0.33 pKTa 1.59 0.65

Comput. LFERb

Computational LFERb

Experimental Computational LFERb a b c d

pK2

pK4

6.80 6.29 ± 0.56

19.14 12.81

pK1,expc

pK2,expd

5.96 [36] 5.20 5.59

13.7 [37] 20.74 13.51

NH2OH M +NH3O. Estimation from linear free energy relationships, Ref. [9]. K1,exp = (K1 + K2). 1 1 K 2;exp ¼ ðK 1 . 3 þ K4 Þ

164

M.I. Fernández et al. / Chemical Physics Letters 490 (2010) 159–164

c.a. KT = 101.6, i.e. NH2OH is favoured, contrary to what happens in simple amino acids, where the zwitterion is the major species by far (c.a. 104) [30]. This result implies the amount of +NH3O in aqueous solution is large enough to play a relevant role in kinetic processes, which is also in agreement with recent structure–activity studies [9]. 4. Conclusions The estimation of several microscopic acid–base equilibria, in aqueous solution, involving NH2OH has been carried out using thermodynamic cycles. The corresponding Gibbs free energies in the gas phase were computed with the CBS-QB3 method, and solvation was modelled by using the conductor-like polarizable continuum model (CPCM). Results depend heavily on the solute cavity, Pauling radii working best. pKa values involving protonated, zwitterionic and neutral species are in reasonable agreement with experiment, whereas the result for anions is worse. The value obtained for the tautomeric equilibrium between NH2OH and +NH3O, KT = ([+NH3O]/[NH2OH]) = 2.6  102, indicates that although NH2OH is favoured, there is enough +NH3O in aqueous solution to be a kinetically relevant species. Acknowledgements Authors wish to thank Dr. Anthony Kirby for the access to his unpublished material. Thanks are also given to the Xunta de Galicia (Spain) for financial support, and to the Centro de Supercomputación de Galicia (CESGA) for computing facilities. Thanks are also given to reviewers for their helpful suggestions. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.cplett.2010.03.047. References [1] L. Zapała, J. Kalembkiewicz, E. Sitarz-Palczak, Biophys. Chem. 140 (2009) 91. [2] J.R. Pliego Jr., Chem. Phys. Lett. 367 (2003) 145.

[3] A.J. Kirby, D.W. Tondo, M. Medeiros, B.S. Souza, J.P. Priebe, M.F. Lima, F. Nome, J. Am. Chem. Soc. 131 (2009) 2023. [4] A.J. Kirby, B.S. Souza, M. Medeiros, J.P. Priebe, A.M. Manfredi, F. Nome, Chem. Commun. (2008) 4428. [5] Y.S. Simanenko, A.F. Popov, T.M. Prokop’eva, E.A. Karpichev, V.A. Savelova, I.P. Suprun, C.A. Bunton, Russ. J. Org. Chem. 38 (2002) 1286. [6] A.J. Kirby, J.E. Davies, A.S. Brandão, P.F. da Silva, W.R. Rocha, F. Nome, J. Am. Chem. Soc. 128 (2006) 12374. [7] F. Haber, Ber. Dtsch. Chem. Ges. 29 (1896) 2444. [8] M. Bronstrup, D. Schroder, I. Kretzschmar, C.A. Schalley, H. Schwarz, Eur. J. Inorg. Chem. 10 (1998) 1529. [9] A.J. Kirby et al., Chem. Commun. 46 (2010) 1302. [10] A.C. Hopkinson, I.G. Csizmadia, Theoret. Chim. Acta (Berl.) 34 (1974) 93. [11] B.T. Hart, Aust. J. Chem. 29 (1976) 231. [12] B.T. Hart, Aust. J. Chem. 29 (1976) 241. [13] L. Radom, J.S. Binkley, J.A. Pople, Aust. J. Chem. 30 (1977) 699. [14] M.D. Burtzoff, L. Peter, P.A. Lepse, D.Y. Zhang, J. Mol. Struc. (Theochem) 619 (2002) 229. [15] E.L. Øiestad, E. Uggerud, Int. J. Mass Spectrom. 185–187 (1999) 231. [16] F. Angelelli, M. Aschi, F. Cacace, F. Pepi, G. de Petris, J. Phys. Chem. 99 (1995) 6551. [17] D.A. Dixon, J.S. Francisco, Y. Alexeev, J. Phys. Chem. A110 (2006) 185. [18] M.J. Frisch et al., GAUSSIAN 03, Revision C.02, Gaussian Inc., Wallingford CT, 2004. [19] M.R. Nyden, G.A. Petersson, J. Chem. Phys. 75 (1981) 1843. [20] J.A. Montgomery Jr., M.J. Frisch, J.W. Ochterski, G.A. Petersson, J. Chem. Phys. 112 (2000) 6532. [21] J.A. Montgomery Jr., M.J. Frisch, J.W. Ochterski, G.A. Petersson, J. Chem. Phys. 110 (1999) 2822. [22] M. Cossi, N. Rega, G. Scalmani, V. Barone, J. Comput. Chem. 24 (2003) 669. and references therein. [23] F. Weinhold, J.E. Carpenter, in: R. Naaman, Z. Vager (Eds.), The Structure of Small Molecules and Ions, Plenum, 1988, p. 227. and references therein. [24] P.A. Guigere, I.D. Liu, Can. J. Chem. 30 (1952) 948. [25] S. Tsunekawa, J. Phys. Soc. Jpn. 33 (1972) 167. [26] R.D. Bach, A.L. Owensby, C. González, H.B. Schlegel, J.J.W. McDouall, J. Am. Chem. Soc. 113 (1991) 6001. [27] A. Johansson, P.A. Kollman, J.F. Liebman, S. Rothenberg, J. Am. Chem. Soc. 96 (1974) 3750. [28] N. Sadlej-Sosnowska, Theor. Chem. Acc. 118 (2007) 281. [29] D.M. Camaioni, C.A. Schwerdtfeger, J. Phys. Chem. A 109 (2005) 10795. [30] R. Stewart, The Proton Applications to Organic Chemistry, Academic Press, Orlando, 1985. [31] M.-L. Hung, M.L. McKee, D.M. Stanbury, Inorg. Chem. 33 (1994) 5108. [32] P. Boulet, F. Gilardoni, J. Weber, H. Chermette, Y. Ellinger, Chem. Phys. 244 (1999) 163. [33] J.E. Del Bene, M.J. Frisch, K. Raghavachari, J.A. Pople, J. Phys. Chem. 86 (1982) 1529. [34] J. Tyrrell, W. Lewis-Bevan, D. Kristiansen, J. Phys. Chem. 97 (1993) 12768. [35] L. Largo, V.M. Rayón, C. Barrientos, A. Largo, P. Redondo, Chem. Phys. Lett. 476 (2009) 174. [36] P. Lume, P. Lahermo, J. Tummavuori, Acta Chem. Scand. 19 (1965) 2175. [37] M.N. Hughes, H.G. Nicklin, K. Shrimanker, J. Chem. Soc. A (1971) 3485.