Theoretical Prediction of p K a in Methanol: Testing SM8 and SMD

Article pubs.acs.org/JPCB

Theoretical Prediction of pKa in Methanol: Testing SM8 and SMD Models for Carboxylic Acids, Phenols, and Amines Elizabeth L. M. Miguel, Poliana L. Silva, and Josefredo R. Pliego* Departamento de Ciências Naturais, Universidade Federal de São João del-Rei, 36301-160, São João del-Rei, MG, Brazil S Supporting Information *

ABSTRACT: Methanol is a widely used solvent for chemical reactions and has solvation properties similar to those of water. However, the performance of continuum solvation models in this solvent has not been tested yet. In this report, we have investigated the performance of the SM8 and SMD models for pKa prediction of 26 carboxylic acids, 24 phenols, and 23 amines in methanol. The gas phase contribution was included at the X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level. Using the proton exchange reaction with acetic acid, phenol, and ammonia as reference species leads to RMS error in the range of 1.4 to 3.6 pKa units. This finding suggests that the performance of the continuum models for methanol is similar to that found for aqueous solvent. Application of simple empirical correction through a linear equation leads to accurate pKa prediction, with uncertainty less than 0.8 units with the SM8 method. Testing with the less expensive PBE1PBE/6-311+G** method results in a slight improvement in the results.

1. INTRODUCTION Organic molecules can be activated for reactions in solution phase through ionization. Among the possible ionization reactions, protonations and deprotonations have a preeminent importance and it is a ubiquitous process in solution chemistry and biochemistry.1−3 This distinguished reaction is associated with the equilibrium constant Ka or its minus log, the pKa. The role of this property for understanding and predicting the solution chemistry has induced a lot of research on theoretical pKa calculations in aqueous solution.4−20 In the present, it is possible to predict pKa of organic molecules with error less than 1 pKa unit provided that a carefully chosen calculation process is used. However, while aqueous solution has received a lot of attention, organic solvents are less studied.21−27 In particular, we have not found any theoretical study of pKa in methanol. Many organic molecules have low solubility in aqueous solution, but are readily soluble in methanol, which is a polar and protic solvent. Methanol has many similarities to water. For example, its ability to solvate cations and anions is close to that of water. Thus, the solvation free energy of the pair of ions NH4+ + Cl− is −158.9 kcal mol−1 in water28 and −156.6 kcal mol−1 in methanol.29 As a consequence, methanol is also able to promote ionization reactions of neutral species like water does.30 On the other hand, the methyl group facilitates the solubilization of less polar organic molecules in methanol. These properties make methanol useful for promoting reactions involving neutral molecules and highly polar transition states. In the present time, methanol is a widely used solvent for chemical reactions. The development of solvation models remain a very active research area.31−48 It is important to have available models that are reliable and easily applied for predicting solvation effects. © 2014 American Chemical Society

Continuum models meet the requirement of practicality and the SMx family of Cramer and Truhlar is parametrized for more than 100 solvents.32,33,49−51 Although these models have been scrutinized for several problems and have presented an overall good performance,52−54 more extensive tests are needed in order to reach the limits of the model. In particular, continuum models are less reliable for solvents with high specific interactions, like hydrogen bonds.15,55 In this way, solvent effects in methanol are an important test for solvation models. The aim of the present study is evaluate the ability of the SM8 and SMD models to predict solvent effects for proton exchange reactions in methanol. The test involves the calculation of pKa for carboxylic acids, phenols, and amines using isodesmic reactions. This part of the study shows how reliable the model is for reactions in solution involving similar species. Other part of the study is related to accurate pKa prediction. Thus, we have also investigated empirical procedures that could be added to the theoretical calculation of pKa for predicting pKa values with uncertainty around 0.5 pKa units in this important solvent.

2. EQUATIONS FOR CALCULATION OF PKA There are different equations for theoretical calculation of pKa. The most direct approach would be the ionization reaction, forming the proton as indicated below: HA → H+ + A−

(1)

Received: February 7, 2014 Revised: May 10, 2014 Published: May 12, 2014 5730

dx.doi.org/10.1021/jp501379p | J. Phys. Chem. B 2014, 118, 5730−5739

The Journal of Physical Chemistry B

Article

Figure 1. Gas-phase free energy of the proton exchange reaction calculated for (a) 20 carboxylic acids with acetate ion, (b) 15 phenols with phenoxide ion, and (c) 11 protonated amines with ammonia. All results at X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level. Experimental data taken from NIST database.76

gas phase free energy contribution, and the second term is the solvation contribution. In the case of cationic acids, the exchange reaction is

Nevertheless, this approach requires the solvation free energy of the proton, a value with considerable uncertainty and long dispute in the literature for aqueous solution.56−66 In methanol, the uncertainty is even greater.29,60 Furthermore, the solvation free energy has the bulk contribution and the surface potential contribution. Thus, solvation models can adopt different solvation free energy scales for calibration of their parameters. In the case of SM8 and SMD models, the solvation free energies of ions and even zwitterionic species are very different. For example, the solvation free energy calculated for the NH4+ ion is −89.7 kcal mol−1 at SM8 level and −79.3 kcal mol−1 at SMD level (calculated in this work). Another approach is the proton exchange scheme, indicated in the equation below: HA + Ref − → A− + HRef

BH+ + Ref → B + RefH+

and the respective equation for pKa calculation is pK a(BH+) =

* ΔGsol + pK a(HRef) RT ln(10)

(2)

pK a(corrected) = a + b*pK a(X3LYP)

(6)

(7)

and the values of a and b are determined by fitting theoretical pKa(X3LYP) values with experimental pKa values. The slope of this correlation is b and the intercept is a. The ideal situation is b close to one and a close to zero, which means that theoretical pKa(X3LYP) values are accurate. Equations 3, 6, and 7 were used in this work to predict more reliable pKa values.

(3)

and * = ΔGg* + ΔΔGsolv * ΔGsol

* ΔGsol + pK a(RefH+) RT ln(10)

The calculation of pKa through eqs 3 and 6 is an interesting test of the solvation model. A reliable model should produce a straight line with slope one and intercept zero. However, due to the uncertainty in the calculated solvation free energy, a simple linear regression can be used to correlate the theoretical and experimental pKa values. Such a procedure could produce even more accurate pKa values. In this way, the corrected theoretical pKa would be calculated by

where HRef is the reference Bronsted acid. Considering that the difference in the solvation free energy of ions should have less uncertainty, this approach should lead to a small error. In addition, this approach is a test of the solvation model. In fact, if the solvation model is able to evenly treat the different species, it must result in accurate pKa. The respective equations are pK a(HA) =

(5)

3. QUANTUM CHEMISTRY CALCULATIONS Molecular structures of the neutral species were initially generated using the MM2 force field in the CHEM3D program. The obtained structures were fully optimized at

(4)

The first term in the left side of eq 4 is the solution phase free energy of reaction in eq 2, the first term on the right side is the 5731



Article

Figure 2. Theoretical versus experimental pKa values for carboxylic acids. Theoretical values at X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level. Data available in Table 1 and in the Supporting Information.

density functional theory level, using the X3LYP functional67,68 and the Dunning DZV+P(d) basis set,69 followed by harmonic frequency calculations. Polarized double-ζ valence basis sets are enough to provide good geometries for minimum energy structures. Analysis of different conformations was done to choose the most stable structure. More accurate electronic energy was calculated using the X3LYP functional and the reliable TZVPP basis set70 extended with diffuse sp functions on the heavy atoms but carbon. In the solvation free energy calculations, we have used the B3LYP/6-31G(d) level for the SM8 method.51 This model is parametrized for specific levels of theory and we have used the parameters for the B3LYP/631G(d) level. For the SMD solvation model,50 the calculations were done with the X3LYP functional and the 6-31(+)G(d) basis set. This basis set is equivalent to the 6-31G(d) basis set for hydrogen and carbon, and the 6-31+G(d) basis set for the remaining atoms. We have chosen the X3LYP functional because it is stable and shows good performance for proton affinity calculations.67,71 The investigated compounds are carboxylic acids (26 species), phenols (24 species), and amines (23 species) with experimental pKa values taken from Rived et al.72 Additional single point energy calculations were done at the PBE1PBE/6-311+G** level using the X3LYP/DZV+P(d) geometries. These additional calculations are based on the findings of Toomsalu et al.,71 which suggest that the PBE1PBE/ 6-311+G** level of theory predicts reliable acidities and basicities when compared with a set of many different functionals. All the geometry optimizations and frequency and single point energy calculations were done with the Firefly program.73

For the solvation, the SM8 calculations were done with the GAMESSPLUS program74 and the SMD calculations were done with a recent version (2011) of the GAMESS program.75

4. RESULTS AND DISCUSSION The present study aims to evaluate the ability of the SM8 and SMD models to predict pKa in methanol solution. Therefore, it is important to determine the reliability of the two terms in eq 4, the gas phase (ΔGg) and the solvation contribution (ΔΔGsolv) to the pKa. In Figure 1, we have provided the values of ΔGg for proton exchange reactions 2 and 5 calculated at X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level and compared with the experimental data taken from NIST database.76 For carboxylic acids, we can notice a good correlation between theory and experiment, with RMS error of only 1.87 kcal mol−1. We have calculated a similar RMS error for phenols, only 1.91 kcal mol−1. This error translates to an uncertainty of 1.4 units in the theoretical pKa value. However, we should notice that many of the experimental gas phase data points have an uncertainty of 2 kcal mol−1. Considering that we are calculating the free energy of isodesmic reactions, it is possible that error cancelation produces more reliable theoretical values. For amines, the performance of theory is better, with an RMS error of only 1.47 kcal mol−1, leading to uncertainty in the pKa of only 1.1 units. Overall, the present findings indicate the X3LYP functional is reliable for computation of the gas phase contribution (ΔGg) to pKa. In addition, our results for the gas phase proton exchange reaction are in agreement with the findings of Toomsalu et al.71 in their evaluation of different functionals for calculating acidities and basicities. Those authors have found better 5732



Article

Table 1. Theoretical and Experimental pKa Values of Carboxylic Acids in Methanola SM8 acetic acid 2-bromoacetic acid 2-chloroacetic acid 2-fluoroacetic acid 2-cyanoacetic acid 2-phenylacetic acid 2,2-dichloroacetic acid propanoic acid 2,3-dichloropropanoic acid benzoic acid 2-fluorobenzoic acid 2-bromobenzoic acid 2-chlorobenzoic acid 2-nitrobenzoic acid 3-chlorobenzoic acid 3-nitrobenzoic acid 3-(trifluoromethyl)benzoic acid 3-cyanobenzoic acid 4-fluorobenzoic acid 4-bromobenzoic acid 4-chlorobenzoic acid 4-nitrobenzoic acid 4-methylbenzoic acid 4-cyanobenzoic acid 2,4-dichlorobenzoic acid 2,4-dinitrobenzoic acid 2,6-dinitrobenzoic acid

SMD

SM8

pKa (X3LYP)

pKa (corrected)

pKa (X3LYP)

pKa (corrected)

pKa (PBE1PBE)

pKa (corrected)

Exp.

8.33 7.29 6.38 6.32 11.90 3.61 11.08 5.77 10.42 9.68 10.37 9.20 7.86 9.62 9.41 9.24 9.83 10.81 10.08 10.25 9.20 11.10 9.75 9.21 5.30 4.61

8.07 7.59 7.17 7.14 9.72 5.90 9.34 6.89 9.03 8.69 9.01 8.47 7.85 8.66 8.57 8.49 8.76 9.21 8.88 8.96 8.47 9.34 8.73 8.48 6.67 6.36

7.12 6.51 6.00 4.75 10.19 2.63 10.81 4.47 9.35 8.25 8.10 7.07 6.71 8.46 7.58 8.14 7.84 9.49 8.94 9.06 7.54 10.11 7.88 6.93 2.99 3.49

8.17 7.91 7.68 7.14 9.51 6.22 9.77 7.02 9.14 8.66 8.59 8.15 7.99 8.75 8.37 8.61 8.48 9.20 8.96 9.01 8.35 9.47 8.50 8.09 6.38 6.60

8.15 7.05 6.17 6.09 11.36 3.36 10.96 5.61 9.69 9.02 9.96 8.76 7.25 8.89 8.77 8.50 9.16 10.04 9.44 9.58 8.64 10.43 9.10 8.78 4.94 3.97

8.23 7.70 7.28 7.24 9.76 5.94 9.57 7.01 8.96 8.64 9.09 8.52 7.79 8.58 8.52 8.39 8.71 9.13 8.84 8.91 8.46 9.32 8.68 8.53 6.69 6.23

9.63 8.06 7.88 7.99 7.50 9.43 6.38 9.71 7.50 9.30 8.41 8.19 8.31 7.64 8.83 8.32 8.69 8.53 9.23 8.93 9.09 8.34 9.51 8.42 7.80 6.45 6.30

Values at 25 °C. X3LYP corresponds to calculations at X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level. PBE1PBE corresponds to PBE1PBE/6311+G**//X3LYP/DZV+P(d) level. The SMD calculations have used the X3LYP/6-31(+)G(d) density and the SM8 calculations have used the B3LYP/6-31G(d) method. The corrected values for each method are explained in the text. Experimental data was taken from ref 72. More details in Supporting Information.

a

performance of the X3LYP functional for basicities than for acidities. The calculated pKa of 26 carboxylic acids versus the experimental values are presented in Figure 2 and Table 1. More detailed data are presented in the Supporting Information. The pKa values range 3.4 units. Acetic acid (pKa = 9.63) was used as the reference species for the proton exchange reaction. Using the same functional group, we should expect a small error in the pKa. However, as we can notice, there is a substantial deviation for both SM8 and SMD methods. We should draw the attention for the trends in the deviations. Stronger acids (lower pKa) have negative deviations, whereas weaker acids (higher pKa) have positive deviations. This is an indication that the continuum solvation models used in this report are not able to treat evenly the different species. Stronger acids have more charge dispersion on the anion form and the solvation is overestimated. The inverse occurs for weaker acids. The RMS errors for this calculation are 1.44 and 1.58 pKa units for SM8 and SMD methods, respectively. The highest deviations are 2.77 and 3.75 pKa units for 2,2dichloroacetic acid at SM8 and SMD levels. Details of the results are presented in the Supporting Information. Although the calculation of pKa with proton exchange reaction using a similar species as the reference molecule has produced an error higher than expected, it is possible to correct the results through a simple linear fit. In the case of carboxylic acids, the correction equations are

pK a(corrected) = 4.230 + (0.4610)· pK a(X3LYP) (SM8, 26 carboxylic acids, RMS = 0.39)

(8)

pK a(corrected) = 5.080 + (0.4341)· pK a(X3LYP) (SMD, 26 carboxylic acids, RMS = 0.21)

(9)

and the respective values of pKa are presented in Figure 2b and d, as well as in Table 1. As we can see, the corrected theoretical values are very close to the experimental data, reaching an RMS error of only 0.21 pKa units for the SMD model. In the case of the SM8 method, the error is also small, 0.39 pKa units. These results are very meaningful, indicating that pKa of carboxylic acids can be predicted accurately in methanol solution using the present approach. Considering the trend in the gas phase ΔGg, the solution phase calculated pKa indicates deficiency of the continuum models to evenly treat the solvation of the different anions. These results point out that the continuum models overestimate the solvation of anions with high charge dispersion (less solvated) and underestimate the solvation of more solvated anions. Further support for this view is the uncertainty in the experimental gas phase acidities and the solution phase pKa. While the former has an uncertainty of 2 kcal mol−1, experimental pKa values have uncertainty of 0.1 unit. The theoretical ΔGg values were calculated for a proton exchange reaction using the same functional group and there is a clear 5733



Article

Figure 3. Theoretical versus experimental pKa values for phenols. Theoretical values at X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level. Data available in Table 2 and in the Supporting Information file.

level of theory presented in this paper in conjunction with the simple correction equations presented below:

correlation between theoretical and experimental data (Figure 1a). On the other hand, the theoretical pKa values need a linear fit with slope as small as 0.43 to predict correct values. With these observations taken together, we can claim that the main source of uncertainty is the ΔΔG*solv term. Phenols were the next class of compounds investigated (Figure 3, Table 2, 24 compounds) and the molecules were selected to have a pKa range of 13 units. We have chosen the parent phenol (pKa = 14.33) as the anchor species. The calculated ΔGg values have a RMS error similar to carboxylic acids, 1.91 kcal mol−1. In fact, there is a clear correlation between theory and experiment as we can see in Figure 1b. On the other hand, in the calculation of pKa, the deviations are relatively high. For the SM8 model, the RMS error is 1.77 pKa units (2,6-dinitrophenol not included due to convergence problem in SM8), whereas for SMD the RMS error reaches 2.56 pKa units. The highest deviation is observed for the 2,4,6trinitrophenol, which amount to 4.40 and 5.71 pKa units for SM8 and SMD methods, respectively. In free energy terms, the error reaches a value as high as 7.8 kcal mol−1. This value is unacceptably high and points out serious deficiency of the continuum solvation models in methanol. The behavior of the theoretical values is similar to carboxylic acids, i.e., more acid species present more negative deviations. On the other hand, when a simple linear fit is used (eqs 10 and 11), the error decreases substantially (Figure 3b and d). The RMS errors become 0.63 and 0.68 pKa units for SM8 and SMD methods, respectively, and the highest deviations decreases to 1.79 pKa units (4-nitrophenol) at both SM8 and SMD levels. Therefore, it is possible to predict reliable pKa values for phenols using the

pK a(corrected) = 4.658 + 0.6658· pK a(X3LYP) (SM8, 23 phenols, RMS = 0.63)

(10)

pK a(corrected) = 5.691 + 0.6025· pK a(X3LYP) (SMD, 24 phenols, RMS = 0.68)

(11)

The third group investigated was the amines, with a pKa range of 12 units. In this case, the correlation between theory and experiment for gas phase basicity is even better (Figure 1c) and the RMS error is 1.47 kcal mol−1. However, for liquid phase pKa values, the deviation increases substantially (Figure 4 and Table 3). Thus, in the case of the SM8 method, the RMS error is 3.63 and decreases to 2.78 at SMD level. Somewhat unexpected, the performance inverts for the corrected values. In fact, although the SMD model has a better performance than SM8 for theoretical calculation of pKa, when the empirical correction is applied, the SM8 model presents a RMS error of only 0.78, while it reaches 1.35 at SMD level. Deviations of the SM8 calculations are more uniform and a simple correction leads to substantial improvement. On the other hand, the SMD model is more erratic and not able to predict accurate pKa values as found for phenols and carboxylic acids. The maximum deviation found for the corrected SM8 method was the value of pyridine, 1.52 pKa units. The equations used for correction are pK a(corrected) = 4.128 + 0.6402· pK a(X3LYP) (SM8, 23 amines, RMS = 0.78) 5734

(12)



Article

Table 2. Theoretical and Experimental pKa Values of Phenols in Methanola SM8 phenol 1-naphtol 2-chlorophenol 2-fluorophenol salicylaldehyde 2-nitrophenol 2-methylphenol 2-methoxyphenol 2-tert-butylphenol 3-bromophenol 3-chlorophenol 3-nitrophenol 3-methylphenol 4-bromophenol 4-chlorophenol 4-nitrophenol 4-methylphenol 4-tert-butylphenol 2,4-dinitrophenol 2,4-dimethylphenol 2-chloro-4-phenylphenol 2,5-dinitrophenol 2,6-dinitrophenol 2,4,6-trinitrophenol 2,4,6-trimethylphenol

SMD

SM8

pKa (X3LYP)

pKa (corrected)

pKa (X3LYP)

pKa (corrected)

pKa (PBE1PBE)

pKa (corrected)

Exp.

12.98 12.57 12.75 11.22 10.23 14.92 16.71 16.13 12.94 12.84 12.02 15.11 13.58 13.61 7.29 15.55 15.68 4.80 16.39 13.05 5.58 −0.85 16.22

13.30 13.03 13.15 12.13 11.47 14.59 15.78 15.40 13.28 13.21 12.66 14.72 13.70 13.72 9.51 15.01 15.10 7.85 15.57 13.35 8.38 4.09 15.46

12.41 11.29 11.73 9.91 9.06 15.05 16.69 16.06 12.35 12.36 10.82 15.00 13.42 13.42 6.33 15.49 15.53 3.33 16.81 11.72 4.19 3.31 −2.16 16.92

13.17 12.49 12.76 11.66 11.15 14.76 15.74 15.37 13.13 13.14 12.21 14.73 13.78 13.77 9.51 15.02 15.05 7.70 15.82 12.75 8.21 7.69 4.39 15.88

12.61 12.65 12.83 11.48 10.51 15.00 17.00 16.20 13.17 12.95 12.20 15.24 13.88 13.93 7.80 15.76 15.92 5.16 16.69 13.20 5.84 −0.76 16.62

12.96 12.98 13.10 12.24 11.62 14.47 15.74 15.24 13.31 13.17 12.69 14.62 13.76 13.79 9.89 14.96 15.06 8.22 15.55 13.33 8.65 4.45 15.50

14.33 13.91 12.97 12.94 12.82 11.53 14.86 14.48 16.50 13.30 13.10 12.41 14.43 13.63 13.59 11.30 14.54 14.52 7.83 15.04 12.70 8.94 7.64 3.55 15.53


a

our RMS error of 1.6 pKa units for 26 carboxylic acids in methanol, although MAD values are smaller than RMS error values. As early as 2002, Klicic et al.10 have applied the Poisson−Boltzmann version of the dielectric continuum model and used direct calculation of pKa in conjunction with linear fit. This approach is mathematically equivalent to proton exchange scheme with linear fit used in this work. They have applied the method for many different functional groups and for carboxylic acids, the MAD was 0.5, close to our RMS error value of 0.4 for SM8 method and less accurate than our RMS error value of 0.2 for SMD. For amines and pyridines, grouped in different sets, the MAD values reported by Sastre et al.18 were 0.4 and 0.8. In our study, these groups were treated altogether and the RMS error was higher, 3.6 and 2.8 pKa units using the SM8 and SMD methods, respectively. Thus, separation into subgroups seems desirable in this case. On the other hand, using the linear fit resolves the problem, leading to RMS error of 0.8 and 1.4 pKa units at SM8 and SMD levels, respectively. Linear fit correction has also been used by Eckert and Klamt13 in the prediction of pKa for organic bases in water using the COSMO-RS solvation method. Those authors have also observed that secondary and tertiary amines require systematic correction of the calculated pKa. When the model including linear fit and systematic correction was applied to a complex and multifunctional set of 58 amines, the RMS error was 0.7 pKa units. This is close to our value for amines using only linear fit and SM8. Sastre et al.18 have also investigated pKa of phenols in aqueous solution using the isodesmic reaction method. The

pK a(corrected) = 2.206 + 0.6446· pK a(X3LYP) (SMD, 23 amines, RMS = 1.35)

(13)

The present results show that the protocol used in this paper can be used to obtain accurate pKa values, with an uncertainty as low as 0.2 for carboxylic acids, 0.6 for phenols, and 0.8 for amines. Methanol and water have similar solvation properties and it is interesting to make a comparison to literature reports. Thus, calculation of pKa through eq 3 using a unique reference species can lead to large error. For example, Pliego and Riveros8 have reported pKa calculations for 15 acids and protonated bases in aqueous solution and considering the hydroxide ion and ammonia as the reference species, respectively. The PCM and SM5.42R methods have presented an RMS error of 7.3 and 7.2 pKa units, respectively. Takano and Houk12 have taken a subset of 11 species from Pliego and Riveros report and the calculations have led to RMS error of 3.0 pKa units with CPCM/UAKS method. When only 7 anions are considered, the RMS error increases to 3.7 pKa units. These results point out how difficult it is for pure continuum solvation models to treat ion solvation. Recently, Sastre et al.18 used the isodesmic proton exchange reaction for pKa calculation in water. The procedure is equivalent to our calculation in this report without linear fitting. They have separated aliphatic and aromatic carboxylic acids and the mean absolute deviation (MAD) was around 1.5 pKa units using B3LYP/6-311++G(d,p) level and solvent effect using the SMD method. This finding is in close agreement with 5735



Article

Figure 4. Theoretical versus experimental pKa values for protonated amines. Theoretical values at X3LYP/TZVPP+diff//X3LYP/DZV+P(d) level. Data available in Table 3 and in the Supporting Information file.

the lower computational cost of this approach. The results are on Tables 1, 2, and 3, and the correction equations are

reported MAD was 0.9 for a set of 10 phenols with a pKa range of only 2.3 units. Our similar calculations for 24 phenols in methanol, involving a pKa range of 13 units, have an error of 1.8 and 2.6 pKa units at SM8 and SMD levels, respectively. With linear fit the error decreases to 0.6 and 0.7 pKa units, respectively. Changing organic solvents from protic to dipolar aprotic can lead to substantial improvement of the continuum solvation models for describing anions.22,23 In the case of cations, Li et al. have reported theoretical calculation of pKa for 98 protonated amines and phosphines in acetonitrile solution.24 Using the PBEPBE/6-311++G(2df,2p)//PBEPBE/6-31G(d) method and PCM solvation, the RMS error was only 1.1 pKa units with slope equal to one. Therefore, theoretical studies of pKa in water, methanol, and dipolar aprotic solvents indicate that protic solvents are a challenging medium for continuum solvation models. The present study makes use of the X3LYP functional with the extended TZVPP+diff basis set. This functional has a better performance than the popular B3LYP method for different properties,67 although for proton affinity the performance is close.71 Even considering that the performance of functionals depends on the corresponding basis set, the use of the triple-ζ basis set is enough to obtain results relatively close to the converged density functional method.70 Thus, the X3LYP/ TZVPP+diff calculations should indicate how accurate the X3LYP functional is. However, the recent report of Toomsalu et al.71 suggesting the good performance of the PBE1PBE/6311+G** method has prompted us to make a test of this level of theory using the SM8 method for solvation. The advantage is

pK a(corrected) = 4.332 + (0.4779) ·pK a(PBE1PBE) (SM8, 26 carboxylic acids, RMS = 0.35)

(14)

pK a(corrected) = 4.939 + (0.6355) ·pK a(PBE1PBE) (SM8, 23 phenols, RMS = 0.62)

(15)

pK a(corrected) = 5.197 + (0.6195) ·pK a(PBE1PBE) (SM8, 23 amines, RMS = 0.75)

(16)

Comparing the X3LYP/TZVPP+diff results, we can notice a slight improvement in the predicted pKa values. In this way, the RMS error decrease from 0.39, 0.63, and 0.78 at X3LYP/ TZVPP+diff level to 0.35, 0.62, and 0.75 at PBE1PBE/6311+G** level for carboxylic acids, phenols, and amines, respectively. Therefore, the present results confirm the good performance of the PBE1PBE/6-311+G** method, and it is another possibility for pKa calculations.

5. CONCLUSION Calculations of proton exchange reaction in gas phase, using the X3LYP/TZVPP+diff//X3LYP/DZ+P(d) method for carboxylic acids, phenols, and amines, have a RMS error of only 1.87, 1.91, and 1.47 kcal mol−1, respectively. These values are within the error bars of the experimental data, indicating that the X3LYP functional with the TZVPP+diff basis set is reliable for obtaining the gas phase reaction free energy. These gas phase 5736



Article

Table 3. Theoretical and Experimental pKa values of Protonated Amines in Methanola SM8

SMD

SM8

pKa (X3LYP)

pKa (corrected)

pKa (X3LYP)

pKa (corrected)

pKa (PBE1PBE)

pKa (corrected)

Exp.

10.88 9.99 7.84 11.49 9.38 4.26 3.78 2.63 −2.92 −1.51 2.86 −7.54 3.52 1.71 −8.87 3.54 −2.81 3.60 3.77 4.43 6.78 10.28 3.68

11.09 10.52 9.14 11.49 10.13 6.85 6.55 5.81 2.26 3.16 5.96 −0.70 6.38 5.22 −1.55 6.39 2.33 6.43 6.54 6.96 8.47 10.71 6.48

12.58 13.66 12.66 13.31 13.77 3.79 4.58 4.71 −0.33 1.31 4.18 −3.77 4.80 3.57 −5.09 5.30 1.63 6.11 5.85 9.03 13.01 14.23 9.85

10.32 11.01 10.37 10.78 11.08 4.65 5.16 5.24 1.99 3.05 4.90 −0.23 5.30 4.51 −1.07 5.62 3.26 6.15 5.98 8.03 10.60 11.38 8.56

10.09 8.76 5.83 10.59 7.84 3.06 2.41 1.04 −4.56 −3.29 1.23 −9.34 1.97 0.06 −10.74 1.98 −4.37 1.99 2.10 2.24 4.70 8.98 1.15

11.45 10.62 8.81 11.76 10.05 7.10 6.69 5.84 2.37 3.16 5.96 −0.59 6.42 5.24 −1.46 6.42 2.49 6.43 6.50 6.58 8.11 10.76 5.91

10.78 11.00 11.20 9.80 11.00 10.78 6.29 5.13 6.05 3.46 3.71 5.95 0.20 5.98 4.95 −0.67 6.57 1.55 7.41 6.89 5.44 7.72 11.07 5.16

ammonia N-methylamine N,N-dimethylamine N,N,N-trimethylamine N-ethylamine N,N,N-triethylamine hydroxilamine o-methylhydroxylamine aniline 2-bromoaniline 2-chloroaniline 2-methylaniline 2-nitroaniline 4-benzylaniline 4-chloroaniline 4-chloro-2-nitroaniline 4-methylaniline 4-nitroaniline 4-hydroxyaniline 4-methoxyaniline pyridine 2,4,6-trimethylpyridine piperidine quinoline


a

data were combined with solvation calculations with SM8 and SMD methods to predict pKa in methanol. The calculated pKa values have an RMS error in the range of 1.4 to 3.6 units when the isodesmic proton exchange reaction is used for carboxylic acids, phenols, and amines. 77 Application of empirical correction through a linear equation leads to substantial improvement in the predicted pKa, with RMS error ranging from 0.2 to 1.4 for SMD and ranging from 0.4 to 0.8 for SM8. Therefore, combination of the SM8 method for solvation contribution with the X3LYP/TZVPP+diff//X3LYP/DZ+P(d) method for gas phase free energy and empirical correction can be used for predicting reliable pKa values in methanol solution. Additional testing with the PBE1PBE/6-311+G**//X3LYP/ DZ+P(d) method for gas phase and SM8 method for solvation leads to a slight improvement in the results.

■

ACKNOWLEDGMENTS

■

REFERENCES

The authors thank the agencies CNPq, FAPEMIG, and CAPES for support.

(1) Anslyn, E. V.; Dougherty, D. A. Modern Physical Organic Chemistry; University Science Books: Sausalito, CA, 2006. (2) Silva, C. M.; Silva, P. L.; Pliego, J. R. Prediction of the pH-rate Profile for Dimethyl Sulfide Oxidation by Hydrogen Peroxide: The Role of Elusive H3O2+ Ion. Int. J. Quantum Chem. 2014, 114, 501−507. (3) Yang, C.; Xue, X.-S.; Jin, J.-L.; Li, X.; Cheng, J.-P. Theoretical Study on the Acidities of Chiral Phosphoric Acids in Dimethyl Sulfoxide: Hints for Organocatalysis. J. Org. Chem. 2013, 78, 7076− 7085. (4) da Silva, C. O.; da Silva, E. C.; Nascimento, M. A. C. Ab Initio Calculations of Absolute pKa Values in Aqueous Solution I. Carboxylic Acids. J. Phys. Chem. A 1999, 103, 11194−11199. (5) Schuurmann, G.; Cossi, M.; Barone, V.; Tomasi, J. Prediction of the pKa of Carboxilic Acids Using the Ab Initio Continuum-Solvation Model PCM-UAHF. J. Phys. Chem. A 1998, 102, 6706−6712. (6) Kallies, B.; Mitzner, R. pKa Values of Amines in Water from Quantum Mechanical Calculations Using a Polarized Dielectric Continuum Representation of the Solvent. J. Phys. Chem. B 1997, 101, 2959−2967. (7) Lim, C.; Bashford, D.; Karplus, M. Absolute pKa Calculations with Continuum Dielectric Methods. J. Phys. Chem. 1991, 95, 5610− 5620. (8) Pliego, J. R., Jr; Riveros, J. M. Theoretical Calculation of pKa Using the Cluster-Continuum Model. J. Phys. Chem. A 2002, 106, 7434. (9) Liptak, M. D.; Shields, G. C. Accurate pKa Calculations for Carboxylic Acids Using Complete Basis Set and Gaussian-N Models

ASSOCIATED CONTENT

S Supporting Information *

Table of gas and solution phase free energy data, as well as theoretical and experimental pKa values. This material is available free of charge via the Internet at http://pubs.acs.org.

■

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 5737



Article

Combined with CPCM Continuum Solvation Methods. J. Am. Chem. Soc. 2001, 123, 7314−7319. (10) Klicic, J. J.; Friesner, R. A.; Liu, S.; Guida, W. C. Accurate Prediction of Acidity Constant in Aqueous Solution via Density Functional Theory and Self-Consistent Reaction Field Methods. J. Phys. Chem. A 2002, 106, 1327−1335. (11) Pliego, J. R., Jr Thermodynamic Cycles and the Calculation of pKa. Chem. Phys. Lett. 2003, 367, 145−149. (12) Takano, Y.; Houk, K. N. Benchmarking the Conductor-Like Polarizable Continuum Model (CPCM) for Aqueous Solvation Free Energies of Neutral and Ionic Organic Molecules. J. Chem. Theory Comput. 2005, 1, 70−77. (13) Eckert, F.; Klamt, A. Accurate Prediction of Basicity in Aqueous Solution with COSMO-RS. J. Comput. Chem. 2006, 27, 11−19. (14) Eckert, F.; Diedenhofen, M.; Klamt, A. Towards a First Principles Prediction of pKa: COSMO-RS and the Cluster-Continuum Approach. Mol. Phys. 2009, 108, 229−241. (15) Ho, J.; Coote, M. L. A Universal Approach for Continuum Solvent pK(a) Calculations: Are We There Yet? Theor. Chem. Acc. 2010, 125, 3−21. (16) Ho, J.; Coote, M. L. pKa Calculation of Some Biologically Important Carbon Acids - An Assessment of Contemporary Theoretical Procedures. J. Chem. Theory Comput. 2009, 5, 295−306. (17) Keith, J. A.; Carter, E. A. Quantum Chemical Benchmarking, Validation, and Prediction of Acidity Constants for Substituted Pyridinium Ions and Pyridinyl Radicals. J. Chem. Theory Comput. 2012, 8, 3187−3206. (18) Sastre, S.; Casasnovas, R.; Muñoz, F.; Frau, J. Isodesmic Reaction for pKa Calculations of Common Organic Molecules. Theor. Chem. Acc. 2012, 132, 1−8. (19) Zhang, S. A Reliable and Efficient First Principles-Based Method For Predicting pKa Values. 4. Organic Bases. J. Comput. Chem. 2012, 33, 2469−2482. (20) Ç iftcioğlu, G. A.; Trindle, C. Computational Estimates of Thermochemistry and pKa Values of Cyclopropenyl Imine Superbases. Int. J. Quantum Chem. 2013, 114, 392−399. (21) Chipman, D. M. Computation of pKa from Dielectric Continuum Theory. J. Phys. Chem. A 2002, 106, 7413−7422. (22) Almerindo, G. I.; Tondo, D. W.; Pliego, J. R., Jr Ionization of Organic Acids in Dimethyl sulfoxide Solution: A Theoretical ab initio Calculation of the pKa Using a New Parametrization of the Polarizable Continuum Model. J. Phys. Chem. A 2004, 108, 166−171. (23) Fu, Y.; Liu, L.; Li, R.; Liu, R.; Guo, Q. First-Principle Predictions of Absolute pKa’s of Organic Acids in Dimethyl Sulfoxide Solution. J. Am. Chem. Soc. 2004, 126, 814−822. (24) Li, J.-N.; Fu, Y.; Liu, L.; Guo, Q.-X. First-Principle Predictions of Basicity of Organic Amines and Phosphines in Acetonitrile. Tetrahedron 2006, 62, 11801−11813. (25) Ding, F.; Smith, J. M.; Wang, H. First-Principles Calculation of pKa Values for Organic Acids in Nonaqueous Solution. J. Org. Chem. 2009, 74, 2679−2691. (26) Eckert, F.; Leito, I.; Kaljurand, I.; Kütt, A.; Klamt, A.; Diedenhofen, M. Prediction of Acidity in Acetonitrile Solution with COSMO-RS. J. Comput. Chem. 2009, 30, 799−810. (27) Trummal, A.; Rummel, A.; Lippmaa, E.; Burk, P.; Koppel, I. A. IEF-PCM Calculations of Absolute pKa for Substituted Phenols in Dimethyl Sulfoxide and Acetonitrile Solutions. J. Phys. Chem. A 2009, 113, 6206−6212. (28) Pliego, J. R., Jr; Riveros, J. M. Gibbs Energy of Solvation of Organic Ions in Aqueous and Dimethyl Sulfoxide Solutions. Phys. Chem. Chem. Phys. 2002, 4, 1622−1627. (29) Pliego, J. R.; Miguel, E. L. M. Absolute Single-Ion Solvation Free Energy Scale in Methanol Determined by the Lithium ClusterContinuum Approach. J. Phys. Chem. B 2013, 117, 5129−5135. (30) da Silva, P. L.; Guimaraes, L.; Pliego, J. R. Revisiting the Mechanism of Neutral Hydrolysis of Esters: Water Autoionization Mechanisms with Acid or Base Initiation Pathways. J. Phys. Chem. B 2013, 117, 6487−6497.

(31) Reinisch, J.; Klamt, A.; Diedenhofen, M. Prediction of Free Energies of Hydration with COSMO-RS on the SAMPL3 Data Set. J. Comput.-Aided Mol. Des. 2012, 26, 669−673. (32) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Generalized Born Solvation Model SM12. J. Chem. Theory Comput. 2012, 9, 609−620. (33) Bernales, V. S.; Marenich, A. V.; Contreras, R.; Cramer, C. J.; Truhlar, D. G. Quantum Mechanical Continuum Solvation Models for Ionic Liquids. J. Phys. Chem. B 2012, 116, 9122−9129. (34) Pomogaeva, A.; Thompson, D. W.; Chipman, D. M. Modeling Short-Range Contributions to Hydration Energies with Minimal Parameterization. Chem. Phys. Lett. 2011, 511, 161−165. (35) Knight, J. L.; Brooks, C. L. Surveying Implicit Solvent Models for Estimating Small Molecule Absolute Hydration Free Energies. J. Comput. Chem. 2011, 32, 2909−2923. (36) Fennell, C.; Dill, K. Physical Modeling of Aqueous Solvation. J. Stat. Phys. 2011, 145, 209−226. (37) Pliego, J. R., Jr. Shells Theory of Solvation and the Long-Range Born Correction. Theor. Chem. Acc. 2011, 128, 275−283. (38) de Lima, G. F.; Duarte, H. A.; Pliego, J. R., Jr. Dynamical Discrete/Continuum Linear Response Shells Theory of Solvation: Convergence Test for NH4+ and OH− Ions in Water Solution Using DFT and DFTB Methods. J. Phys. Chem. B 2010, 114, 15941−15947. (39) Liu, J.; Kelly, C. P.; Goren, A. C.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G.; Zhan, C.-G. Free Energies of Solvation with Surface, Volume, and Local Electrostatic Effects and Atomic Surface Tensions to Represent the First Solvation Shell. J. Chem. Theory Comput. 2010, 6, 1109−1117. (40) Mennucci, B. Continuum Solvation Models: What Else Can We Learn from them? J. Phys. Chem. Lett. 2010, 1, 1666−1674. (41) Tan, C.; Tan, Y. H.; Luo, R. Implicit Nonpolar Solvent Models. J. Phys. Chem. B 2007, 111, 12263−12274. (42) Boes, E. S.; de Andrade, J.; Stassen, H.; Goncalves, P. F. B. Computational Study of Anion Solvation in Nitrobenzene. Chem. Phys. Lett. 2007, 436, 362−367. (43) Dzubiella, J.; Swanson, J. M.; McCammon, J. A. Coupling Nonpolar and Polar Solvation Free Energies in Implicit Solvent Models. J. Chem. Phys. 2006, 124, 084905. (44) Soteras, I.; Curutchet, C.; Bidon-Chanal, A.; Orozco, M.; Luque, F. J. Extension of the MST Model to the IEF Formalism: HF and B3LYP Parametrizations. THEOCHEM 2005, 727, 29−40. (45) Nerenberg, P. S.; Jo, B.; So, C.; Tripathy, A.; Head-Gordon, T. Optimizing Solute−Water van der Waals Interactions to Reproduce Solvation Free Energies. J. Phys. Chem. B 2012, 116, 4524−4534. (46) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3093. (47) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Dill, K. A. Predictions of Hydration Free Energies from All-Atom Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 4533−4537. (48) Nicholls, A.; Wlodek, S.; Grant, J. A. The SAMP1 Solvation Challenge: Further Lessons Regarding the Pitfalls of Parametrization. J. Phys. Chem. B 2009, 113, 4521−4532. (49) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on the Generalized Born Approximation with Asymmetric Descreening. J. Chem. Theory Comput. 2009, 5, 2447− 2464. (50) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (51) Marenich, A. V.; Olson, R. M.; Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Self-Consistent Reaction Field Model for Aqueous and Nonaqueous Solutions Based on Accurate Polarized Partial Charges. J. Chem. Theory Comput. 2007, 3, 2011−2033. (52) Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Prediction of SAMPL2 Aqueous Solvation Free Energies and Tautomeric Ratios Using the SM8, SM8AD, and SMD Solvation Models. J. Comput.-Aided Mol. Des. 2010, 24, 317−333. 5738



Article

(73) Granovsky, A. A. Firefly v 7.1.F; http://classic.chem.msu.su/ gran/gamess/index.html. (74) Higashi, M.; Marenich, A. V.; Olson, R. M.; Chamberlin, A. C.; Pu, J.; Kelly, C. P.; Thompson, J. D.; Xidos, J. D.; Li, J.; Zhu, T.; et al. GAMESSPLUS, v 2010−2; University of Minnesota, Minneapolis, MN, 2010. (75) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (76) Linstrom, P. J.; Mallard, W. G. NIST Chemistry WebBook; NIST Standard Reference Database Number 69, 2011. (77) For amines, some reactions are not truely isodesmic.

(53) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation Modeling. Acc. Chem. Res. 2008, 41, 760−768. (54) Chamberlin, A. C.; Cramer, C. J.; Truhlar, D. G. Performance of SM8 on a Test to Predict Small-Molecule Solvation Free Energies. J. Phys. Chem. B 2008, 112, 8651−8655. (55) Pliego, J. R., Jr; Riveros, J. M. The Cluster-Continuum Model for the Calculation of the Solvation Free Energy of Ionic Species. J. Phys. Chem. A 2001, 105, 7241−7247. (56) Vlcek, L.; Chialvo, A. A.; Simonson, J. M. Correspondence Between Cluster-Ion and Bulk Solution Thermodynamic Properties: on the Validity of the Cluster-Pair-Based Approximation. J. Phys. Chem. A 2013, 117, 11328−11338. (57) Beck, T. L. The Influence of Water Interfacial Potentials on Ion Hydration in Bulk Water and Near Interfaces. Chem. Phys. Lett. 2013, 561−562, 1−13. (58) Donald, W. A.; Williams, E. R. An Improved Cluster Pair Correlation Method for Obtaining the Absolute Proton Hydration Energy and Enthalpy Evaluated with an Expanded Data Set. J. Phys. Chem. B 2010, 114, 13189−13200. (59) Fawcett, W. R. The Ionic Work Function and its Role in Estimating Absolute Electrode Potentials. Langmuir 2008, 24, 9868− 9875. (60) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Single-Ion Solvation Free Energies and the Normal Hydrogen Electrode Potential in Methanol, Acetonitrile, and Dimethyl Sulfoxide. J. Phys. Chem. B 2007, 111, 408−422. (61) Lamoureux, G.; Roux, B. Absolute Hydration Free Energy Scale for Alkali and Halide Ions Established from Simulations with a Polarizable Force Field. J. Phys. Chem. B 2006, 110, 3308−3322. (62) Grossfield, A.; Ren, P. Y.; Ponder, J. W. Ion Solvation Thermodynamics from Simulation with a Polarizable Force Field. J. Am. Chem. Soc. 2003, 125, 15671−15682. (63) Asthagiri, D.; Pratt, L. R.; Ashbaugh, H. S. Absolute Hydration Free Energies of Ions, Ion-Water Clusters, and Quasichemical Theory. J. Chem. Phys. 2003, 119, 2702−2708. (64) Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R. The Proton’s Absolute Enthalpy and Gibbs Free Energy of Solvation From ClusterIon Solvation Data. J. Phys. Chem. A 1998, 102, 7787−7794. (65) Pliego, J. R., Jr; Riveros, J. M. On the Calculation of the Absolute Solvation Free Energy of Ionic Species: Application of the Extrapolation Method to Hydroxide Ion in Aqueous Solution. J. Phys. Chem. B 2000, 104, 5155−5160. (66) Hunenberger, P. H.; Reif, M. M. Computation of MethodologyIndependent Single-Ion Solvation Properties from Molecular Simulations. III. Correction Terms for the Solvation Free Energies, Enthalpies, Entropies, Heat Capacities, Volumes, Compressibilities, and Expansivities of Solvated Ions. J. Chem. Phys. 2011, 134, 144103. (67) Xu, X.; Zhang, Q.; Muller, R. P.; Goddard, W. A., III An Extended Hybrid Density Functional (X3LYP) with Improved Descriptions of Nonbond Interactions and Thermodynamic Properties of Molecular Systems. J. Chem. Phys. 2005, 122, 014105−014114. (68) Xu, X.; Goddard, W. A., III The X3LYP Extended Density Functional for Accurate Descriptions of Nonbond Interactions, Spin States, and Thermochemical Properties. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 2673−2677. (69) Dunning Jr, T. H.; Hay, P. J. In Methods of Electronic Structure Theory; Schaefer, H. F., III, Ed.; Plenum Press: New York, 1976; pp 1−27. (70) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (71) Toomsalu, E.; Koppel, I. A.; Burk, P. Critical Test of Some Computational Chemistry Methods for Prediction of Gas-Phase Acidities and Basicities. J. Chem. Theory Comput. 2013, 9, 3947−3958. (72) Rived, F.; Rosés, M.; Bosch, E. Dissociation Constants of Neutral and Charged Acids in Methyl Alcohol. The Acid Strength Resolution. Anal. Chim. Acta 1998, 374, 309−324. 5739


Theoretical Prediction of p K a in Methanol: Testing SM8 and SMD

Recommend Documents