1. Introduction

The important role played a part by proton transfer reactions in chemistry, bioorganic chemistry, biology, and material science[1,2] makes the assessment of reliable pKa values, a topic of continuing connotation. The deprotonation energies of organic compounds and the proton affinities of the corresponding conjugate bases are widely used for the prediction of gas-phase and aqueous phase Bronsted acidities[3–10]. Strong acids have small values of deprotonotion energy while strong bases have large values of proton affinity. Several works on the prediction of the acidity of organic compounds can be found in the literature. For instance, Smith and Radom[11,12] have shown that the G2 and G2(MP2) methods provide excellent results for both deprotonation enthalpies and proton affinities of small molecules and Liton et al. have successfully calculated pKa values of trimethylaminium ion using G2(MP2) gas phase Gibb’s free energy together with B3LYP/6-311G(d)-CPCM solvation energy[13]. Catalan and Palomar[14] have investigated gas phase acidities of a number of species and have shown that calculations at the B3LYP method with 6-311+G(d) and 6-311+G(3df,3pd) basis sets correlate well with the experimental data. Good correlations have been obtained between experimental pKa values of a wide range of organic compounds and their calculated gas-phase deprotonation enthalpies[15].

Correlations between theoretical predictions and gas-phase acid-base equilibrium constants of organic compounds have been reported for amines, alcohols and thiols[16–18]. An excellent correlation was obtained between the aqueous-phase acidity calculated with the HF/3-21G(d) method and experimental pKa values for a series of nitrogen bases[19]. Recently, various theoretical descriptors were used to investigate their correlation with the carbon acidity of compounds. A reasonable correlation was obtained between the deprotonation energies of the compounds calculated at the HF/3-21G(d) and B3LYP/6-31G(d) levels of theory and their aqueous pKa values[20,21]. The ability to predict acidity using a coherent, well-defined theoretical approach, without external approximation or fitting to experimental data would be very useful to chemists. However, the current situation is less satisfactory in solution, mostly due to the difficulty of calculating solvation energies with adequate accuracy. The aim of this work is to establish the theoretical calculation without any external approximation or fitting for determining the pKa values of some organic compounds that exhibit pKa values from 5 to 50.

2. Method and Theoretical Calculations

2.1. Thermodynamic Cycle

The gas-phase Gibb’s free energy change ∆Gº_gas and aqueous phase Gibb’s free energy change ∆Gº_aq in Figure 1 are calculated using Eq. (2) and Eq. (4)[23] respectively. The equilibrium constant of reaction (1) is Ka and the pKa is -logKa for this reaction.

Table 1. The following compounds have been studied for theoretical pKa calculations

There is a relation between total change of Gibb’s free energy in Figure 1, ∆G°_aq and pKa. Acidity constant (pKa) was calculated using the thermodynamic cycle[24] shown in Figure 1. In this cycle, ΔG°_gas and ΔG°_aq are respectively, the gas-phase and aqueous-phase standard Gibb’s free energy change. ΔG°(s, HA), ΔG°(s, A^-) and ΔG°(s, H⁺) are also the standard free energies of solvation for HA, A^- and H⁺ respectively.

(1)

(2)

(3)

(4)

(5)

(6)

(7)

Computationally calculating pKa constants are a demanding and laborious process because an errors of 1.36 kcal/mol in the change of free energy in Eq. (1) results in an error of 1 pKa unit. There are at least three sources of error in pKa calculations[25]. The first is the model used to calculate pKa, which generally involves a thermodynamic cycle shown in Figure 1 and the second and third major errors stem from the calculations of ∆Gº_gas and ∆G°(s). In this work, gas-phase Gibb’s free energies and the change in free energies of solvation (aqueous) of HA and A^- have been calculated accurately to minimize errors.

Figure 1. Thermodynamic cycle used to calculate ∆Gº aq using Eq. (4) from its components: ∆Gºgas is the Gibb’s free change in the gas-phase; ∆Gº(s, HA), ∆Gº (s, A-) and ∆Gº(s, H+) are the solvation energies of HA, A- and H+ respectively at 298.15ºK

2.2. Computational Details

All the calculations were performed on Intel core-i5 Fujitsu Laptop Computer with 3-4 GB of memory and 48 GB of scratch disk space, using Gaussian98 Revision-A.9[26] software. Every structure is initially optimized in the gas phase at HF/6-31G(d,p) level of theory. The geometries of the neutral and its anions are finally optimized at the B3LYP/6-311+G(2d,2p) and MPW1PW91/6-311+ G(2d,2p) methods and frequency calculations are performed at the same level of theory to compute the ∆G°_gas values (Gas phase Gibb’s free energy change) as well as the stationary points were characterized with no imaginary frequency. We also used HF and B3LYP methods at 6-31G(d,p) and 6-311++G(2d,2p) basis sets respectively with the conductor-like polarizable continuum model (CPCM)[27] for the calculation of ∆G°(s) (Gibb’s free energy change of solvation). Before the calculation of solvation energies all the structure were finally optimized in the presence of solvent (water) at HF/6-31+G(d) level of theory and the optimized structures were confirmed to be energy minima by vibrational frequency calculations with all real frequencies. These levels of theory have been designed particularly for the prediction of reliable energies of molecules in the gas-phase and have been demonstrated to provide the accuracy within 1–2 kcal/mol when assessed against large test sets of thermo-chemical data[28, 29].

Table 2. Gibb’s free energy of each molecule and its anion together with Gibb’s energy changes of reaction (1), ∆Gºgas, in the gas phase at 298.15ºC

3. Results and Discussion

3.1. Selection and Validation of Experimental Data

First we re-examined some reported experimental and theoretical values, which we use in this work for comparison or as part of the calculations. The solvation energy of H⁺ in water is very important in the calculation of pKa values. However, to date, several theoretical and experimental values have been reported: -264.0[30], -264.6[31], -265.9[32, 33] and -266.1[32] kcal/mol. This resulting pKa value can then be compared with the (directly-measured) experimental pKa value of the compounds. The value of -265.9 kcal/mol was supported by Tissandier et al.[34], on the basis of experimental ion water clustering data and the aqueous solvation free energies of neutral ion pairs, after the correction of standard states for the gas-phase and aqueous solution. Thus, they reported a value of -264.0[30] kcal/mol for the absolute aqueous solvation free energy of the proton and by the correction of 1.89 kcal/mol[34] for converting the standard state of 1 atm to 1 mol/L, the value of -265.9 kcal/mol can be obtained[32]. In this article the recent experimental gas phase Gibb’s free energy of -6.28 kcal/mol[22] and the change in free energy of solvation (aqueous) of -264.0 kcal/mol[30] of H⁺ at standard state are used to calculate the pKa values.

3.2. Calculation of the pKa

The procedure sketched above has been applied to the calculation of the pKa values for a set of nineteen (19) organic compounds shown in Table 1. The compounds 1-7 are designated as carbon compounds or carbon acids and the rest 8-19 are alcoholic compounds or oxygen acids. Some of these molecules have been successfully studied by Mekelleche and his group[23]. In this paper numerous methods have been attempted to compute the accurate gas phase Gibb’s free energy changes for pKa values but the most effective methods are B3LYP/6-311+G(2d,2p) and MPW1PW91/6-311+G(2d,2p) level of theory. Gibb’s free energy of each molecule and its anion in the gas phase together with Gibb’s free energy changes of reaction (1), ∆Gº_gas and the experimental values with unsigned errors are shown in Table 2. It is observed that as oxygen is directly attached to the compounds, the calculation of gas phase Gibb’s free energy favours MPW1PW91/6-311+G(2d,2p) level of theory. It may be due to the electrostatic effects arising from unbalanced delocalized charges. Comparison between calculated energies and experimental values as well as unsigned errors show that the calculated ∆Gº_gas values of carbon acids (1-7) at B3LYP/6-311+G(2d,2p) and the values for oxygen acids (8-19) at MPW1PW91/6-311+G(2d,2p) level of theory are very good agreement to the experimental values.

Different types of methods have been employed in order to calculate accurate solvation energies. Comparatively the most successful methods are HF/6-31G(d,p) and B3LYP/6-311++G(2d,2p) level of theory for the prediction of solvation energies. Now both methods are employed to compute the pKa values shown in Table 3 and 4. The pKa values obtained from the solvation energy at B3LYP/6-311++G(2d,2p)-CPCM method combined with the gas phase energy at B3LYP/6-311+G(2d,2p) method for carbon acids (1-7) are very close to the experimental values shown in Table 3. On the other hand, the pKa calculated from the solvation energy at HF/6-31G(d,p)-CPCM method combined with the gas phase Gibb’s free energy at MPW1PW91/6-311+G(2d,2p) method for oxygen acids (8-19) are also very close to the experimental values shown in Table 4. Most of the calculated solvation energies are very close to the experimental values for neutral molecules but for anion species a small deviation exists. As expected, both models provide comparable results for neutral molecules, whereas larger differences are obtained for anionic species, where escaped charge effects become more significant. It may be due to the strongly solvated in solution and the strong electrostatic effects arising from unbalanced delocalized charges. Inspection of calculated pKa reveals that in most cases the values are reasonably very close to the experimental values. For example, the calculated pKa are 32.04, 16.07 and 14.40 for the compounds 2, 10 and 15 whereas the experimental pKa are 32.0, 16.10 and 14.42 respectively.

Table 3. Calculated gas and aqueous-phase ∆Go values of the acid dissociation together with solvation energy (in kcal/mol) as well as calculated and experimental aqueous pKa values at 298.15 ºK

In other cases, though small deviations are observed but it is not unforeseen. Figure 2 shows the correlation between the experimental and calculated pKa values with lower unsigned errors shown in Table 3 and 4. The lower unsigned errors indicate the accuracy of calculated pKa values. A useful way to consider agreement between experiment and theory is to calculate the linear regression between the theoretical and experimental quantities. The degree of agreement between the two is then reflected by the slope and intercept of the correlation line, the standard deviation from the correlation line and the correlation coefficient. The last two depend on the scatter of the points about the correlation line. A slope different from unity and a non-zero intercept imply systematic deviation between experiment and theory and means that the theoretical values of relative acidities will be systematically either too high or too low. A non-zero mean difference between theory and experiment reflects an overall bias in the absolute values of the calculated acidity[7]. The regression equation with the corresponding standard deviation (SD) and square correlation coefficient (R²) values are given in Eq. (8). The F- and P-value are 30890.46 and 0.000, respectively, which confirm the validity of regression results.

Figure 2. Correlation of the experimental pKa with the calculated pKa values

(8)

From the correlation equation and Figure 2 we observe that the slope is very near to unity and the line passes through the origin, the predicted pKa are in close proximity to the experimental values.

Table 4. Calculated gas and aqueous-phase ∆Go values of the acid dissociation together with solvation energy (in kcal/mol) as well as calculated and experimental aqueous pKa values at 298.15 ºK

4. Conclusions

In this paper we have employed HF and DFT (B3LYP, MPW1PW91) methods combined with conductor like polarized continuum dielectric (CPCM) aqueous solvation calculations to designate acidities of 19 compounds. We have successfully calculated accurate pKa for some carbon acids and oxygen acids. The pKa calculated from the gas phase energy at B3LYP/6-311+G(2d,2p) combined with the solvation energy at B3LYP/6-311++G(2d,2p)-CPCM methods for the carbon acids and also the gas phase energy at MPW1PW91/6-311+G(2d,2p) combined with solvation energy at HF/6-31G(d,p)-CPCM methods for oxygen acids show an excellent agreement to the experimental pKa values. The values obtained from these methods could sometimes approximately match to the experimental values while in some cases small deviations are found. Finally, we have attempted to develop an effective method for determining pKa with the highest degree of agreement with experimental values using low cost methods with small basis sets instead of high level of theory as well as without any external approximation or fitting.