Theoretical calculation of pKas of phosphoric (V) acid in the polarisable continuum and cluster-continuum models

Journal of Molecular Structure 924–926 (2009) 170–174

Contents lists available at ScienceDirect

Journal of Molecular Structure journal homepage: www.elsevier.com/locate/molstruc

Theoretical calculation of pKas of phosphoric (V) acid in the polarisable continuum and cluster-continuum models Maciej S´miechowski * ´ sk University of Technology, Narutowicza 11/12, 80-952 Gdan ´ sk, Poland Department of Physical Chemistry, Chemical Faculty, Gdan

a r t i c l e

i n f o

Article history: Received 1 September 2008 Received in revised form 5 November 2008 Accepted 12 November 2008 Available online 14 December 2008 Keywords: Phosphoric acid pKa Polarisable continuum model Cluster-continuum model Ab initio calculations

a b s t r a c t The accuracy of the polarisable continuum model (PCM) and the mixed cluster-continuum model in the prediction of the absolute values of the three consecutive pKas of phosphoric (V) acid has been checked. PCM calculations at the MP2/6-31+G(d,p) level reproduce the first pKa with an acceptable error. However, they fail significantly for the next two pKas and increasing the level of theory to G3B3 compound method does not provide any improvement. On the other hand, cluster-continuum calculations at the same MP2 level adequately predict all three dissociation constants of H3PO4. The number of necessary solvating water molecules depends on the polarising power of the anion and increases with the charge of the phosphate group. The obtained results indicate the validity of the cluster-continuum approach for the prediction of reliable pKa values of polyprotic inorganic acids. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction First principles determination of pKa of Brønsted acids and bases in an aqueous solution has received a considerable attention in the recent years [1–5]. The ‘‘natural” state of reactants in the computational methods, the gas phase, needs to be amended, however, to provide a reasonable description of the solutes in the continuous solvent phase. This is usually formalised as an interaction of the solute embedded in a chosen cavity with the surrounding dielectric medium. The definition of the cavity and the details of its construction lead to a number of different approaches to the problem. These range from the simple spherical cavity (Onsager) model [6,7], through a cavity defined as a sum of fixed overlapping spheres, as in the polarisable continuum model (PCM) [8–10], finally to a ‘‘molecular shape” cavity defined by an electronic isodensity surface (IPCM and SCIPCM models) [11]. The simplest spherical cavity approximation predicts zero solvation energy for systems lacking dipole moment (like the PO3 4 anion encountered in this study), which limits its applicability. The PCM approach, especially in the new parameterisation proposed for the purpose of calculation of accurate solvation free energies (the so-called UAHF atomic radii set) [12], is able to predict ionisation constants of simple monoprotic acids in aqueous solution with an acceptable error [3]. However, it usually fails to provide a reliable description of protonation equilibria for more complex systems, like polyprotic acids, very strong or very weak * Tel.: +48 58 347 1283; fax: +48 58 347 2694. E-mail address: [email protected] 0022-2860/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.molstruc.2008.11.047

acids, or zwitterionic molecules [1,4,5]. The situation might be significantly improved by explicit consideration of several solvent molecules around the ionisable group of the solute. This removes the most important drawback of the dielectric continuum approach – the lack of directionality of solute–solvent interactions. This problem is especially important in the highly-structured solvents, e.g., water. A formal description of the solvation of such aqueous clusters by a continuum has been provided by Pliego and Riveros in their cluster-continuum model [1,13]. To date, it was successfully applied in calculations of solvation free energies [13], pKas of simple acids [1], ionisation constants of biological zwitterionic molecules [4,5] and in the investigation of reaction pathways in aqueous solutions [14,15]. In the original formulation, the discrete water-solute cluster is embedded in an isodensity-based cavity, as defined by the IPCM model [1,11]. The phosphate anions and their organic derivatives are crucially important in biological systems [16]. Particularly, the H2 PO 4 and anions play a key role in the metabolic pathways, since at HPO2 4 physiological pH these are the most abundant forms present in an aqueous solution [17]. The phosphate buffer system, comprised of both these anions, is one of the major pH controlling systems in the human body, particularly important in the intracellular fluid and urine [18]. Ionic equilibria in the phosphoric (V) acid system are complicated and involve multiple species and complexes [17]. This acid, like most inorganic oxoacids, is characterized by dissociation constants rising stepwise by ca. 5 pKa units [19]. Although their values are known experimentally with high certainty, due to their importance, they have to our best knowledge

M. S´miechowski / Journal of Molecular Structure 924–926 (2009) 170–174

never been reproduced by theoretical methods. Such an investigation might be viewed as a trial case for various solvation models used to calculate pKas, since a polyprotic acid with small highly charged anions provides a difficult test, especially for calculating the free energies in solution. In this work, we assessed the effectiveness of the pure continuum model (namely, PCM) vs. the cluster-continuum model in the prediction of three dissociation constants of H3PO4.

The prototypic reaction (1) and its corresponding DGaq define the pKa of an aqueous Brønsted acid. Please note that in this work the HA and A notation does not necessarily imply the charge of  2 the reactants, since HA ¼ H3 PO4 ; H2 PO 4 or HPO4 and A ¼ H2  2 3 PO4 ; HPO4 or PO4 :

HAaq ! Hþaq þ Aaq

ð1Þ

Hereafter DG* denotes Gibbs free energy referred to a standard state of 1 mol L1 as opposed to DG° referring to a standard state of 1 atm. The two standard states for the chemical potential of a species are related by: 

lgas ¼ lgas þ RT lnð0:082053TÞ

ð2Þ

The free energy of reaction (1) in solution is obtained by:

DGaq ¼ DGgas þ DDGsolv

ð3Þ

DGgas

Here is the gas phase reaction free energy and the last term is given by Eq. (4).

¼

DGsolv ðHþ Þ

þ

DGsolv ðA Þ



DGsolv ðHAÞ

ð4Þ

DGsolv (X)

The solvation free energy is calculated according to the Ben–Naim definition [20], so it corresponds to the Xgas(1 mol L1) ? Xaq(1 mol L1) process. Finally, the pKa of the HA acid is obtained from the well-known relation:

pK a ¼ DGaq =ð2:303RTÞ

ð5Þ

Although Eq. (1) provides a simple and elegant definition of pKa, the practical application of this equation is limited by its drawbacks. Firstly, it creates a charge separation from a neutral molecule, an uneasy task for computational chemistry. More importantly, it is based on a free proton – a species nonexistent in any coordinating solvent. The latter limitation might be potentially overcome by considering an explicitly solvated molecule instead (like H3 Oþ aq ), but considering the well-known issue of proton identity in an aqueous solution [21], such an approach seems somewhat arbitrary. Many authors apply the experimental DGsolv (H+) and DGf ðHþgas Þ instead in their calculations [3–5], although the accuracy of values of these parameters raises some doubts [22]. Another alternative is to consider proton transfer from an acid molecule to a hydroxide anion as below:

HAaq þ OHaq ! H2 Oaq þ Aaq

ð6Þ

DGaq

Its is calculated just like for Eq. (1) from Eq. (3) and equation analogous to Eq. (4), but the final pKa calculation must take into account the water dissociation equilibrium that needs to be added to reaction (6) to make it formally equivalent to reaction (1). This leads to Eq. (7) for the pKa of the HA acid.

pK a ¼ DGaq =ð2:303RTÞ þ 15:74

[3], their deficiencies with calculating pKas are also well known and might be potentially overcome by considering discrete aqueous clusters instead of bare reactants [1,4,5]. The extension to the cluster-continuum model herein is based on reaction (6) that, taking into account the hydrating water, might be rewritten as:

HAðH2 OÞm þ OH ðH2 OÞn ! ðm þ n þ 1  pÞH2 O þ A ðH2 OÞp

ð8Þ

Free energy of reaction (8) in solution is calculated from Eq. (3), as before. The correct standard state gas phase and solvation contributions are obtained from Eq. (9) and (10), respectively.

2. Theoretical calculation of pKa

DDGsolv

171

ð7Þ

The procedure outlined above, based on either of the Eqs. (1) or (6), allows the calculation of pKa in the purely continuum models, e.g., PCM. Although they perform quite satisfactorily in many cases



DGgas ¼ DGgas þ ðm þ n  pÞRT lnð0:082053TÞ DDGsolv

¼

ð9Þ

ðm þ n þ 1  pÞDGsolv ½H2 O þ DGsolv ½A ðH2 OÞp   DGsolv ½HAðH2 OÞm   DGsolv ½OH ðH2 OÞn 

ð10Þ

Comparing reactions (6) and (8), we see that the equilibrium constant of the latter differs from that of the former only by the H2O concentration. Taking that into account, we finally arrive at the proper formula for calculating pKa of the HA acid in the cluster-continuum approach:

pK a ¼ DGaq =ð2:303RTÞ þ 15:74 þ ðm þ n  pÞ log½H2 O

ð11Þ 1

with the H2O concentration at 298.15 K equal to 55.344 mol L . In this work, we use the continuum model based on reaction (1) or (6) with pKa given by Eq. (5) or (7), respectively. In turn, the cluster-continuum approach utilizes Eqs. (8)–(11), and (3) to obtain the pKa values. 3. Determination of the optimal cluster size in the clustercontinuum model 

Let us define the ‘‘incremental water binding energy” ðDEn Þ as the energy of the process of adding a water molecule to an aqueous cluster of a solute S.

SðH2 OÞn1 þ H2 O ! SðH2 OÞn

ð12Þ

Of course, this reaction can be considered both in the gas phase and in an aqueous solution, giving DEn;gas and DEn;aq , respectively (note that DEn ¼ DEn , since the choice of the ideal gas standard state is irrelevant for DE or DH at a constant temperature). In the context of this work, DEn;aq is more meaningful, as the added water molecule should stabilize the cluster in an aqueous solution. A given water molecule, to be added to a cluster, must favourably interact with it. That is, DEn;aq for its addition should be more negative than the water–water interaction energy. The latter may be approximated by the negative energy of vaporization of H2O, or simply by its energy of solvation. 



DEint ½H2 O ¼ DEvap ½H2 O ¼ DEsolv ½H2 O

ð13Þ

The above considerations form the grounds for the determination of the optimal number of water molecules in an aqueous cluster for the purpose of application in the cluster-continuum model. Namely, the largest n value, for which the condition that DEn;aq < DEsolv [H2O] is still fulfilled, is considered as the optimal number of H2O molecules to be used in Eq. (8). 4. Computational methodology The continuum calculations were performed according to the polarisable continuum model (PCM) [8–10] as implemented in Gaussian 03. Two independent approaches consisted of: (a) MP2/ 6-31+G(d,p) [23–27] optimization and harmonic frequency analysis in the gas phase, followed by a HF/6-31+G(d,p) calculation of DGsolv using the PCM formalism, UAHF atomic radii set [12] and

M. S´miechowski / Journal of Molecular Structure 924–926 (2009) 170–174

172

SCFVAC option of Gaussian 03; (b) G3B3 compound method calculation in the gas phase [28], followed by a HF/GTLarge calculation of DGsolv as above. GTLarge is an exhaustive basis set defined in the last step of the G3B3 method [28]. The cluster-continuum calculations took account of the bulk solvent by applying the isodensity polarisable continuum model (IPCM) [11]. The first step was calculation of the gas phase optimal structures of aqueous clusters and their harmonic vibrational frequencies at the MP2/6-31+G(d,p) level. All structures were verified to be true local energetic minima by the absence of imaginary frequencies. Solvation contribution to free energy was then calculated according to the IPCM formalism with the MP2/6-31+G(d,p) wave function using an isodensity of 0.0004 and a dielectric constant of 78.39, as proposed before [1,13]. Frozen core approximation was not used in the MP2 calculations. The MP2/6-31+G(d,p) level of theory applied here is generally considered to be the smallest model giving adequate structures of hydrogen-bonded species [29]. All zero-point energies obtained at this level were multiplied by a recommended scale factor of 0.9657 [30]. All calculations were performed with the Gaussian 03 system [31]. 5. Results and discussion Although it was argued that aqueous solution geometry is sometimes necessary to predict more reliable pKas, the actual difference in pKa derived using the gas phase or the aqueous solution geometry is mostly <1 unit [3]. Furthermore, the optimization of ionic clusters using PCM with the UAHF atomic radii set often results in unphysical geometry of the system, as checked in the preparation of this work. On the other hand, the IPCM approach allows single point calculations only. Therefore, reoptimization of the clusters in the continuum was eventually not attempted and gas phase geometries only were used for the purpose of calculating DGsolv values. The obtained ab initio results in the PCM formalism are summarized in Tables 1 and 2. Considering the individual reactants, it is immediately seen that, although the solvation free energy of the hydroxide ion is predicted extremely well (cf. experimental DGsolv of 447 kJ mol1 [32]), the same conclusion cannot be claimed for the phosphate anions. The experimental DGsolv values [32] for 1 1 ) and PO3 ) are much more H2 PO 4 (481 kJ mol 4 (2781 kJ mol negative than the calculated ones. The inadequacy of the PCM UAHF model for the prediction of reliable solvation free energies for complex and multiple charged anions was noted before [4]. In

Table 1 Calculated gas-phase and solvation free energies of the reactants in the H3PO4 dissociation pathway in the PCM approximation. Reactant

MP2/6-31+G(d,p)a 

DGgas +e

H H2O OH H3PO4 H2 PO 4 HPO2 4 3 PO4 a

b

0.010012 76.232606 75.612541 642.917669 642.401232 641.679693 640.762170

G3B3a 

DGsolv c

DGgas b

DGsolv d

1112.48 30.67 442.96 61.55 267.69 933.41 2081.71

0.010012 76.401395 75.777956 643.943321 643.422162 642.692321 641.769603

1112.48 27.78 442.37 56.82 262.04 928.81 2089.91

Level of theory for the gas phase calculations. Gas phase free energy (1 atm standard state, hartree). c Solvation free energy at the HF/6-31+G(d,p) level in the PCM UAHF model (1 mol L1 standard state, kJ mol1). d Solvation free energy at the HF/GTLarge level in the PCM UAHF model (1 mol L1 standard state, kJ mol1).  e Experimental values taken from Ref. [33] ðDGgas Þ and Ref. [34] ðDGsolv Þ. b

Table 2 Calculated reaction free energies in aqueous solution and pKa values for the consecutive H3PO4 dissociation steps in the PCM approximation. Step

e

I If IIe IIf IIIe IIIf a b c d e f

MP2/6-31+G(d,p)a

DGaq b

pKa

26.79 65.92 105.70 12.99 137.66 44.95

3.30 4.20 17.13 18.02 22.73 23.62

c

G3B3a

DpKa

d

1.14 2.04 9.92 10.81 11.79 10.41

DGaq b

pKac

DpKad

40.11 59.16 126.45 27.18 138.50 39.23

5.64 5.38 20.76 20.50 22.87 22.61

3.48 3.22 13.55 13.29 10.55 10.29

Level of theory for the gas phase calculations. Reaction free energy in aqueous solution calculated from Eq. (3) (kJ mol1). pKa value calculated from Eq. (5) or (7). Deviation from experimental value (Ref. [19]). Results based on prototypic reaction (1). Results based on prototypic reaction (6).

the PCM calculations we used the experimental results for H+: DGogas ¼ 26:29 kJ mol1 [33] and DGsolv ¼ 1112:48 kJ mol1 [34]. Such approach is often argued to lower the calculated pKa errors [3–5]. The obtained pKa values (Table 2) clearly reveal the failure of the pure continuum models to predict the further dissociation constants of polyprotic acids. The overall performance of reaction (6) and reaction (1) for pKa calculations was similar. Interestingly, the G3B3 compound method was found much inferior to the less demanding MP2 calculations with moderate basis set. Since the solvation contributions for both methods are practically the same, except for the last dissociation step (cf. Table 1, the resulting difference in DDGsolv is lower than 3.5 kJ mol1 for steps I and II), the larger part of the error must come, quite surprisingly, from the gas phase reaction free energies. The free energies of individual reactants are more negative at the G3B3 level, as expected (Table 1), but the resulting reaction free energies are more positive at this level, leading to larger values of DGaq and pKa (Table 2). The experimental values of pKas of H3PO4 are 2.16, 7.21 and 12.32 [19]. The first pKa was found to be reproduced with an acceptable accuracy, especially in the MP2 calculations. The deviation from experiment was within the error limits of the PCM UAHF method [3]. However, pKa,2 and pKa,3 were predicted with an error exceeding 10 units. The conclusion that the polarisable continuum models perform poorly in the more difficult cases was reached before, especially for ions with charge localised on oxygen atoms, as for the phosphate family investigated here [1]. The cluster-continuum (or ‘‘microsolvated”) models require the optimized structures of aqueous clusters to be known. We therefore performed respective calculations for the following numbers of solvating water molecules: OH(H2O)1–4, H3PO4(H2O)1–3, 2 3 H2 PO 4 (H2O)1–4, HPO4 (H2O)1–6 and PO4 (H2O)1–6. As starting points in the optimizations, we used the most energetically favourable structures found before for aqueous clusters of OH [35], 2 3 H2 PO 4 [36], HPO4 [37] and PO4 [38]. The incremental water binding energies DEn;aq necessary for the determination of the optimal number of solvating water molecules for each solute were calculated from the cluster thermal energies in the IPCM approximation and they are shown in Table 3. In theory, they should be compared to the negative experimental energy   of vaporization of water ðDEvap ½H2 O ¼ DHvap ½H2 O RT = 41.51 kJ 1 at 298.15 K [39]). However, more consistent approach mol would be to compare them to the computed solvation energy of water in IPCM, which is equal to 22.90 kJ mol1. Applying the criterion that DEn;aq < DEsolv [H2O], the following optimal numbers of water molecules in each cluster were found: OH(H2O)3, 2 3 H3PO4(H2O)2, H2 PO 4 (H2O)1, HPO4 (H2O)4 and PO4 (H2O)6. For the latter anion, even at n = 6 the requirement is still fulfilled by a wide margin. However, optimization of larger clusters was not at-

M. S´miechowski / Journal of Molecular Structure 924–926 (2009) 170–174 Table 3 Incremental water binding energies of the studied clusters in an aqueous solution.a nb

OH

H3PO4

H2 PO 4

HPO2 4

PO3 4

1 2 3 4 5 6

–67.98 46.70 34.79 18.08

28.97 25.86 11.08

26.64 21.28 18.37 11.17

49.99 47.02 44.99 38.10 19.12 16.09

97.03 79.88 71.49 68.22 58.23 47.33

a b

DEn;aq

as defined by Eq. (12) (kJ mol ). Number of water molecules in the cluster.

H3 PO4 ðH2 OÞ2 þ OH ðH2 OÞ3 ! H2 PO4 ðH2 OÞ þ 5H2 O

ð14Þ

H2 PO4 ðH2 OÞ þ OH ðH2 OÞ3 ! HPO2 4 ðH2 OÞ4 þ H2 O 2  3 HPO4 ðH2 OÞ4 þ OH ðH2 OÞ3 ! PO4 ðH2 OÞ6 þ 2H2 O

ð15Þ ð16Þ

The reagents’ and reaction free energies necessary for calculation of the three pKas are summarized in Tables 4 and 5. The solvation free energies of the bare ions calculated with the IPCM methodology, included for comparison with the values for clusters, are actually much worse than their PCM counterparts (Table 4). Though conceptually very intuitive, the isodensity approach to defining the molecular cavity is known, however, to perform poorly with solvation energies of charged species in an aqueous solution [3]. The strength of the cluster-continuum approach comes from the accurate modelling of the gas-phase as well as solvation contributions to reaction (8), by inclusion of discrete solvent

Table 4 Calculated gas-phase and solvation free energies of the reactants in the H3PO4 dissociation pathway in the cluster-continuum approximation. 



Reactant

na

DGgas b

DGsolv c

DGsolv;n¼0 d

H2O OH H3PO4 H2 PO 4 HPO2 4 3 PO4

– 3 2 1 4 6

76.232606 304.378326 795.391345 718.645527 946.718878 1098.439175

22.90 225.48 55.64 232.40 670.97 1383.79

– 291.15 47.59 246.37 848.76 1789.72

a

Number of water molecules in the cluster. Gas phase free energy of the cluster (1 atm standard state, hartree). Solvation free energy of the cluster in the IPCM model (1 mol L1 standard state, kJ mol1). d Solvation free energy of the bare ion in the IPCM model (1 mol L1 standard state, kJ mol1). b

c



Step

Eq.a

DGaq b

pKac

DpKad

I II III

14 15 16

136.16 45.99 18.98

1.14 7.69 14.16

3.30 0.48 1.84

a b

d

tempted here, as the PO3 4 (H2O)6 represents the whole first hydration sphere of PO3 4 in the gas phase [38]. Pliego and Riveros noted previously that most often the ideal number of solvent molecules varies between two and three [1]. The values found here differ from this observation, but this apparent disagreement might be rationalised taking into account the properties of the hydrated ions. H2 PO 4 is a large monovalent anion with a small polarising power (q/r = 0.5, where q is the charge number and r is the ionic radius in Å [40]), so it is supposed to  be rather weekly hydrated, similarly to, e.g., ClO4 [41]. In turn, 2 3 HPO4 and PO4 , due to their high charge, have large polarising powers (q/r = 1.0 and 1.26, respectively, basing on experimental ionic radii [40]), which implies strong hydration in an aqueous solution, exceeding even the sulphate case [41]. Very large hydration spheres for phosphate anions have been found recently in neutron diffraction experiments [42]. Considering the optimal number of water molecules needed to stabilize the reacting species, Eq. (8) for the three dissociation steps of H3PO4 might be rewritten explicitly in the form of Eqs. (14)–(16) below:



Table 5 Calculated reaction free energies in aqueous solution and pKa values for the consecutive H3PO4 dissociation steps in the cluster-continuum approximation.

c

1

173

Equation number in the text. Reaction free energy in aqueous solution calculated from Eq. (3) (kJ mol1). pKa value calculated from Eq. (11). Deviation from experimental value (Ref. [19]).

molecules ‘‘buffering” the solute from the dielectric continuum. In the highly-structured solvents, such as water, the directionality of intermolecular interactions cannot be neglected. In the clustercontinuum approach, the solvating water molecules improve the description of this ‘‘directional” part of the interactions (comprised almost exclusively of hydrogen bonding), whereas the outside continuum adequately represents the long-range part of the interactions, since liquid water is known to lack long-range ordering and only coulombic forces act at such distances [43]. The pKa values in Table 5 demonstrate the advantage of the cluster-continuum approach over the continuum models without explicit consideration of the solvent. The improvement in the prediction of pKa,2 and pKa,3 is impressive. The model overestimates the first dissociation constant of H3PO4, but the error in its determination is nonetheless comparable to the continuum results. The following pKas are underestimated, as are all the PCM results, but the error with respect to experiment is minimized considerably. The root-mean-square error for pKa determination found previously for a large testing data set was equal to 2.2 units [1]. Our results (Table 5) are but in a single case lower than this benchmark value. The plausible explanation of the success of the cluster-continuum model is the lowering of the DDGsolv contribution to Eq. (3), which is still more difficult to reproduce accurately than the DGgas value [1]. Since most of the solvation energy is ‘‘hidden” within the cluster, the dielectric continuum represents the lesser part of it. This seems particularly important for highly charged anions studied in this work, for which DGsolv increases quickly with increasing polarising power. Let us note finally that the tedious task of optimizing multiple cluster structures is not always necessary, since very often the optimal number of water molecules can be assumed by analogy to the cases studied in more detail [1,13]. 6. Conclusions This study compares the accuracy of the cluster-continuum model as opposed to the simple polarisable continuum model in the prediction of the three consecutive pKas of phosphoric (V) acid. In accordance with previously published results for simple acids [3], we have found that the PCM model with the UAHF parameterisation of ionic radii [12] is adequate for reproducing the experimental value of the first pKa of H3PO4. Ab initio calculations at the nowadays computationally inexpensive MP2/6-31+G(d,p) level are sufficient for this purpose and even appear superior to the G3B3 compound method coupled with PCM/UAHF determination of solvation free energy. However, the continuum methods fail considerably to predict the two following pKas of H3PO4. An almost twofold overestimation of pKa,3 indicates an error of the same order in the prediction of DGaq of prototypic reaction (1), revealing the inaccuracy of the PCM approach for thermochemical calculations in solution in the more difficult cases. The extension of the solvation model to include explicitly the nearest water molecules significantly improves the description of

M. S´miechowski / Journal of Molecular Structure 924–926 (2009) 170–174

174

the phosphoric acid dissociation reactions. All three dissociation constants of H3PO4 could be adequately reproduced with the IPCM variation of the cluster-continuum model [1,13]. The optimal number of solvating water molecules increases with the charge of the phosphate group. The obtained results indicate the validity of the cluster-continuum approach for the prediction of reliable pKa values of polyprotic acids and open this area for future studies. Especially, the extension of this work towards biologically important molecules with phosphate functional groups (e.g., O-phospho-Lserine) seems interesting. Acknowledgments This work was supported from the Republic of Poland scientific funds as a research project, within Grant No. N N204 3799 33. Calculations were carried out at the Academic Computer Centre in Gdan´sk (TASK). The author is indebted to Professor Janusz Stangret for his on-going support and many fruitful discussions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

J.R. Pliego, J.M. Riveros, J. Phys. Chem. A 106 (2002) 7434. J.R. Pliego, Chem. Phys. Lett. 367 (2003) 145. G.A.A. Saracino, R. Improta, V. Barone, Chem. Phys. Lett. 373 (2003) 411. J. Pratuangdejkul, W. Nosoongnoen, G.-A. Guérin, S. Loric, M. Conti, J.-M. Launay, P. Manivet, Chem. Phys. Lett. 420 (2006) 538. M. Canle, D.R. Ramos, J.A. Santaballa, Chem. Phys. Lett. 417 (2006) 28. L. Onsager, J. Am. Chem. Soc. 58 (1936) 1486. M.W. Wong, M.J. Frisch, K.B. Wiberg, J. Am. Chem. Soc. 113 (1991) 4776. S. Miertus, E. Scrocco, J. Tomasi, Chem. Phys. 55 (1981) 117. S. Miertus, J. Tomasi, Chem. Phys. 65 (1982) 239. M. Cossi, V. Barone, R. Cammi, J. Tomasi, Chem. Phys. Lett. 255 (1996) 327. J.B. Foresman, T.A. Keith, K.B. Wiberg, J. Snoonian, M.J. Frisch, J. Phys. Chem. 100 (1996) 16098. V. Barone, M. Cossi, J. Tomasi, J. Chem. Phys. 107 (1997) 3210. J.R. Pliego, J.M. Riveros, J. Phys. Chem. A 105 (2001) 7241. J.R. Pliego, J.M. Riveros, Chem. Eur. J. 8 (2002) 1945. J.R. Pliego, Chem. Phys. 306 (2004) 273. W. Saenger, Principles of Nucleic Acid Structure, Springer-Verlag, New York, 1984 (p. 82).

[17] J. Baril, J.-J. Max, C. Chapados, Can. J. Chem. 78 (2000) 490. [18] B.M. Koeppen, Adv. Physiol. Educ. 275 (1998) 132. [19] D.R. Lide (Ed.), CRC Handbook of Chemistry and Physics, 88th ed, CRC Press, Boca Raton, 2008 (pp. 8-40–8-41). [20] A. Ben-Naim, J. Phys. Chem. 82 (1978) 792. [21] G.A. Voth, Acc. Chem. Res. 39 (2006) 143. [22] J.R. Pliego, J.M. Riveros, Chem. Phys. Lett. 332 (2000) 597. [23] C. Møller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [24] R. Ditchfield, W.J. Hehre, J.A. Pople, J. Chem. Phys. 54 (1971) 724. [25] W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. [26] V.A. Rassolov, J.A. Pople, M.A. Ratner, T.L. Windus, J. Chem. Phys. 109 (1998) 1223. [27] T. Clark, J. Chandrasekhar, G.W. Spitznagel, P.v.R. Schleyer, J. Comput. Chem. 4 (1983) 294. [28] A.G. Baboul, L.A. Curtiss, P.C. Redfern, K. Raghavachari, J. Chem. Phys. 110 (1999) 7650. [29] J.E. Del Bene, M.J.T. Jordan, J. Mol. Struct. (Theochem) 573 (2001) 11. [30] J.P. Merrick, D. Moran, L. Radom, J. Phys. Chem. A 111 (2007) 11683. [31] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery Jr., T. Vreven, K.N. Kudin, J.C. Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford CT, 2004. [32] Y. Marcus, Ion Properties, Marcel Dekker Inc., New York, 1997 (pp. 120–129). [33] I.A. Topol, G.J. Tawa, S.K. Burt, A.A. Rashin, J. Chem. Phys. 111 (1999) 10998. [34] M.D. Tissandier, K.A. Cowen, W.Y. Feng, E. Gundluach, M.H. Cohen, A.D. Earhart, J.V. Coe, T.R. Tuttle, J. Phys. Chem. A 102 (1998) 7787. [35] M. S´miechowski, J. Stangret, J. Phys. Chem. A 111 (2007) 2889. [36] C.C. Pye, M.R. Michels, Can. J. Anal. Sci. Spectrosc. 50 (2005) 70. [37] C.C. Pye, M.R. Michels, Can. J. Anal. Sci. Spectrosc. 49 (2004) 175. [38] C.C. Pye, W.W. Rudolph, J. Phys. Chem. A 107 (2003) 8746. [39] E.F.G. Herington, Pure Appl. Chem. 40 (1974) 391. [40] Y. Marcus, Ion Properties, Marcel Dekker Inc., New York, 1997 (pp. 46–54). [41] J. Stangret, T. Gampe, J. Phys. Chem. A 106 (2002) 5393. [42] P.E. Mason, J.M. Cruickshank, G.W. Neilson, P. Buchanan, Phys. Chem. Chem. Phys. 5 (2003) 4686. [43] T. Head-Gordon, G. Hura, Chem. Rev. 102 (2002) 2651.