A Reliable and Efficient First Principles-Based ... - ACS Publications

Dec 4, 2009 - Quantum Solutions, 2013 Green Acres Road, Suite A, FayetteVille, Arkansas 72703. ReceiVed: July 15, 2009; ReVised Manuscript ReceiVed: ...
2 downloads 0 Views 1MB Size
Download PDF
432

J. Phys. Chem. A 2010, 114, 432–442

A Reliable and Efficient First Principles-Based Method for Predicting pKa Values. 2. Organic Acids Shuming Zhang,† Jon Baker,†,‡ and Peter Pulay*,† Department of Chemistry and Biochemistry, UniVersity of Arkansas, FayetteVille, Arkansas 72701, and Parallel Quantum Solutions, 2013 Green Acres Road, Suite A, FayetteVille, Arkansas 72703 ReceiVed: July 15, 2009; ReVised Manuscript ReceiVed: NoVember 10, 2009

In part 1 of this series, we developed a protocol for the large-scale calculation of pKa values in aqueous solutions from first principles calculations, with the goal of striking a compromise between accuracy and computational efficiency. Following previous workers in the field, pKa values are calculated from a linear regression fit to deprotonation energies: pKa(f) ) Rf(EA- - EHA) + βf, where f denotes a family of functional groups. In this paper, we derive (Rf, βf) values for the acidic functional groups -COOH, -POOH, alcoholic and phenolic -OH, -SH, -NHOH/dNOH, and -NROH, using a data set of 449 experimental pKa values. Several groupings of these functional groups were explored; our final recommended method uses five families (10 empirical parameters). Mean absolute deviations between our fits and experiment are 0.4 pKa units or less for each with a maximum error range of (1.5 pKa units. In certain subgroups, such as monocarboxylic acids, considerably better fits (mean absolute deviation ∼0.20 pKs units) were obtained at the cost of more empirical parameters. Almost 70% of pKa’s calculated by our protocol lie within (0.4 pKa units and over 90% within (0.8 pKa units of the experimental reference value. Our results compare favorably with previous similar models which have greater computational cost. 1. Introduction The pKa value of a molecule determines the relative concentration of its protonated and deprotonated forms at a specific pH. Most drugs contain at least one site that is able to protonate or deprotonate reversibly. Protonation or deprotonation greatly affects the physiological binding and transport properties of a molecule. Consequently, the pKa value of a drug provides important information about its efficacy in pharmacokinetics. Experimentally, pKa values (also called acid dissociation constants or acidity constants) can be determined by a variety of methods, including potentiometric titration,1 infrared (IR) and Raman spectrometry,2 conductometry,3 nuclear magnetic resonance (NMR),4 liquid chromatographic-mass spectrometry (LCMS),5 capillary electrophoresis (CE),6 and thermometric methods such as isothermal titration calorimetry (ITC).7 Unfortunately, these are all rather time and sample-consuming with high-cost (particularly NMR, LC-MS, and CE which require a fairly high level of expertise to interpret the results). Theoretical approaches can be in principle much cheaper, offer considerable time savings, and are much easier to use compared with the experimental techniques. There have been literally thousands of papers published in the past decade or so on the theoretical prediction of pKa values using a whole host of computational methodologies. Some examples for organic acids include: carboxylic acids;8-11 alcohols;9-11 benzoic acids;8,10-12 phenols;8,10,13 hydroxamic acids, barbituric acids, and sulfonamides;10 imides;8 thiols;9 uracils;11,14 and tetrazoles.15,16 To be fair, the results have been mixed, but very acceptable mean absolute deviations between theory and experiment ranging from 0.1 to 0.8 pKa units depending on the different functional groups have been obtained.13,17-19 † ‡

University of Arkansas. Parallel Quantum Solutions.

In paper 1 of this series20 we developed a computational protocol for what we hope will be a fast and reliable method for predicting pKa values. For the acid HA, our method involves a full geometry optimization of both HA and the anion A- using the COSMO solvation model (with water as solvent)21 at the OLYP/3-21G(d) level, followed by single-point OLYP22 energy calculations, again with COSMO, using the 6-311+G** basis set.23 We then fit the energy difference between the two species to experimental values via a simple linear regression: pKa(f) ) Rf(EA- - EHA) + βf, where f denotes a family of functional groups. This approach assumes that entropic and other omitted effects can be absorbed into the empirical fitting parameters. In accordance with this, only the energy of the most stable conformer is used. The error introduced by not using a proper Boltzmann average was estimated20 to be about 0.12 pKa units for a typical case when the deprotonated form has only one stable conformer but the protonated form has two. In paper 120 we derived the computational protocol using a training set of 34 small molecules (carboxylic acids, alcohols and phenols, pyridines, anilines, and amines). In this work we derive the actual fitting equations for organic acids using far larger data sets. We confine ourselves here primarily to acidic O-H and S-H groups, leaving the more demanding N-H group to a later publication. 2. Computational Methods All calculations were carried out with the PQS program package24 running on a small PC cluster (800 MHz PIV Xeon processors) or on Red Diamond, a University-wide cluster of 256 64-bit Xeon processors.25 Structures were constructed, job input prepared, and jobs submitted using PQSMol, a graphical user interface and model builder. The full data set used to derive the linear regression fits is provided as Supporting Information, together with a wide range

10.1021/jp9067087  2010 American Chemical Society Published on Web 12/04/2009

Method for Predicting pKa Values of experimental references. It contains 370 different organic acids: carboxylic acids (both aliphatic and aromatic; 82); dicarboxylic acids (33); phosphonic acids (48); aliphatic alcohols (28); phenols (67); thiols (40); hydroxamic acids (45); and oximes (27). Note that some of these compounds (e.g., all the dicarboxylic and phosphonic acids) have two acidic hydrogens, and where this was the case, provided experimental data were available, both pKa values were used in the fits. As a result of this, the data set actually has 449 different pKa values. One compound, 2-phosphonobenzoic acid, occurs in two individual group data sets (carboxylic and phosphonic acids) as the first proton is lost from the carboxyl group (pKa1) while the second comes from the phosphonic group (pKa2). The hydroxamic acids (functional group -NH-OH) were not straightforward, as a proton can be extracted from either the O-H or the N-H group; often, the respective pKa values were very close. We finally decided to treat the hydroxamic acids collectively irrespective of which group lost the proton, as experimentally this would be difficult to determine and furthermore would complicate any subsequent use of our fitting equations. Note that originally we had 47 hydroxamic acids in our data set, but two (4-hydroxybenzohydroxamic acid and 3-hydroxy-2-naphthohydroxamic acid) were removed, as our calculations indicated that the first pKa was actually the loss of a proton from the ring phenol group and it was not clear exactly what was being measured experimentally26,27 Several of the systems in our data set were N-substituted species for which loss of the O-H proton was the only possibility. Unlike the set of 34 small molecules used in paper 1 to determine the theoretical method, the experimental data for the far greater range of compounds considered in this work is much more limited. For many systems we could find only one experimental determination of the pKa; in other cases two (or more) values were found that were significantly different (and hence at least one of them had to be wrong). Several of these are discussed where appropriate later in the text. Thus the experimental data are nowhere near as well established as in paper 1. (Note that all relevant molecules from paper 1 are included in the data sets used here.) Any of these experimental pKa values that require serious revision in the future will of course affect the final linear regression results we derive for the various functional groups. However, given the large size of our data sets and the limited number of fitting parameters, it would require an extensive set of corrections to have a significant effect. A further major issue that caused us far more problems than we initially anticipated was conformational flexibility. We believe that previous workers have also underestimated the importance of this problem. The entropy-free model we use provides for only one computed energy for the neutral form (the acid, HA) and one for the anion (A-); this of course should be the lowest energy conformer. Identifying this beyond a reasonable doubt was not always straightforward. Finding the global minimum among a number of local minima is still an unsolved problem in general, and (when this was not obvious) we did multiple calculations on the same compound using different starting structures. For the molecules in this set where significant doubts arose, we used the GMMX feature of the PCMODEL program28 to explore the configurational space for low-energy (within a 3.5 kcal/mol window) conformers. The GMMX configuration search method was inspired by the Still group.29 It varies randomly the torsional angles around designated bonds, calculates the energy using the MMX force field,28 and saves the lowest energy conformers. The lowest, or a few

J. Phys. Chem. A, Vol. 114, No. 1, 2010 433 lowest conformers, found within 1000 search steps (for simple systems) or 5000 steps (for more complex systems) were subjected to minimization at the default OLYP/3-21G(d)/ COSMO level (unless the conformation was already calculated) and the lowest energy conformer found at this level was used in the higher-level single point calculations. Molecules where this procedure was applied are marked by GMMX in the Supporting Information. Neither this procedure nor any other polynomial-scaling procedure is guaranteed to locate the global minimum. Consequently, although we have a fairly high degree of confidence, we cannot state with complete certainty that we have found the global minimum (the lowest energy conformer) in all cases. It would have been desirable to evaluate the energies of all GMMX low-energy conformers at the DFT level but, in view of its considerable cost, this has not been done. As with the experimental uncertainty, the discovery of lower energy structures for any of the molecules in our data set could slightly alter the final linear regression equations. This was not an issue in paper 120 as either the compounds used there were small or the conformation of the lowest energy structure was known. One further comment: Initially, we divided the compounds for a given ionizable functional group into a training set, which was used to generate preliminary parameters, and a somewhat smaller test set. In all cases the preliminary training-set parameters gave a good fit for the test set. However, given that we always combined the two sets for our final fits, we do not report the intermediate results, except for some discussion for the carboxylic acids. One reason for dividing the data into a training and a test set is to avoid overfitting, which is mainly useful when the number of fitted parameters is commensurate with the number of data points. This is obviously not the case for our fits. 3. Results and Discussion 3.1. Carboxylic Acids. Carboxylic acids, both aliphatic and aromatic (such as benzoic acids), are distinguished from other organic acids because they are typically much easier to deprotonate. The same remarks apply to the phosphonic acids. Experimental pKa values for the carboxylic acids in our data set range from 0.5 to just over 6.0. (See Table 2). Initially we fitted the aliphatic monocarboxylic acids (49) and the aromatic monocarboxylic acids (33) separately. The resulting linear regression plots are shown in Figure 1, and the corresponding linear equations (for these and all other functional groups discussed in this work) in Table 1. The mean absolute deviation between the fitted and experimental pKa values are 0.19 and 0.22 pKa units, respectively (with correlation coefficients (R2) of 0.95 and 0.89). The two carboxylic acid plots are so similar that they can safely be merged into one, giving an almost unchanged mean absolute deviation of 0.20 pKa units (R2 ) 0.93) over 82 compounds. We used the linear fit derived from the 82 monocarboxylic acids to predict pKa values (both pKa1 and pKa2) for a series of 33 dicarboxylic acids (see Table S1 in the Supporting Information.). Under normal circumstances this would not be advisible; it would clearly be more sensible to include some dicarboxylic acids in the training set if one was planning to predict dicarboxylic acid pKa values. However, given that our ultimate aim is to predict pKa’s for drug candidate molecules, many of which will have multiple sites for removal or addition of a proton, this can be considered as a test of how linear regression equations derived principally from molecules with just one (de)protonation site work in systems with several such sites.

434

J. Phys. Chem. A, Vol. 114, No. 1, 2010

Zhang et al.

Figure 1. Carboxylic acids

TABLE 1: Linear Regression Fits to the Various Functional Groups compounds

no. pKa

slope

intercept

R2

MAD

aliphatic carboxylic acids aromatic carboxylic acids all monocarboxylic acids mono- + dicarboxylic acids phosphonic acids carboxylic + phosphonic acids phenols alcohols phenols + alcohols all hydroxamic acids hydroxamic acids (-NHOH) hydroxamic acids (-NROH) oximes hydroxamic acids (-NH-OH) + oximes thiols phenols + alcohols + thiols

49 33 82 148 90 238 67 28 95 45 33 12 29 62 42 137

0.2476 0.2326 0.2428 0.2251 0.2451 0.2366 0.3312 0.3333 0.3295 0.0582 0.2029 0.1223 0.2280 0.2228 0.3769 0.3350

-68.0102 -63.6828 -66.6235 -61.5300 -67.5525 -64.9392 -89.7135 -90.4470 -89.2317 -8.7682 -51.5091 -28.6533 -58.7143 -57.3054 -103.3833 -90.9108

0.95 0.89 0.93 0.85 0.98 0.94 0.96 0.97 0.98 0.08 0.61 0.60 0.95 0.92 0.95 0.98

0.19 0.22 0.20 0.39 0.28 0.36 0.26 0.40 0.30 0.50 0.29 0.17 0.39 0.33 0.36 0.34

Furthermore, given that our final fits will be derived by combining both the training and test sets, then it is an internal test only. Perhaps not surprisingly the results were not particularly impressive when compared to the original monocarboxylic acid fits: R2 was 0.77 with a mean absolute deviation of 0.61 pKa units. However, if all 66 pKa values for the 33 dicarboxylic acids are fitted via linear regression then the resulting linear equation has R2 ) 0.79 and a mean absolute deviation of 0.59 pKa units. This is the best that can be achieved for the dicarboxylic acids alone and is only marginally better than the results from our derived equation. Our rationale for the poorer performance for dicarboxylic acids is that they exhibit specific solvation effects if the two carboxylic groups are spatially close. A continuum solvation model like COSMO is not able to describe such effects. We took all the monocarboxylic (82) and all the dicarboxylic acid (66) pKa’s and divided them into two sets: a training set

error range -0.66 -1.30 -1.37 -1.47 -1.05 -1.53 -1.49 -1.08 -1.51 -1.20 -1.08 -0.48 -0.71 -1.19 -0.81 -1.36

to to to to to to to to to to to to to to to to

+1.00 +0.31 +0.98 +1.69 +0.93 +1.81 +1.30 +1.31 +1.30 +1.10 +0.48 +0.31 +0.71 +0.81 +1.00 +1.38

with 107 pKa’s and a test set with 41 pKa’s. (See the comments under Table S1 in the Supporting Information.) A linear regression fit to the training set gave R2 ) 0.88 with a mean absolute deviation of 0.36 pKa units; using the linear equations derived from the training set to predict pKa values for the test set gave R2 ) 0.77 with a mean absolute deviation of 0.49 pKa units, while a linear regression fit to the test set gave essentially the same results (R2 ) 0.77, mean absolute deviation ) 0.48); the error ranges were also very similar: -1.24 to +1.52 and -1.29 to +1.51 pKa units, respectively. These two sets of results for the test set, predicted from the training set and fitted to the test set, are extremely close, showing that a good choice of training set can indeed provide very accurate predictive linear regression equations. It also indicates that we need a wide range of compounds, including those with more than one (de)protonation site, if we wish to obtain good predictive linear regression equations for drug candidate molecules.

Method for Predicting pKa Values If we fit the experimental acid dissociation constants of all mono- and dicarboxylic acids (115 compounds; 148 pKa values) with a single line, the resulting fit has R2 ) 0.85 and a mean absolute deviation of 0.39 pKa units (Figure 1d). This is of course very much in line with the partial fits discussed above and represents our final predictive linear regression equation for carboxylic -OH groups in drug candidate molecules. If the predicted pKa for a simple monocarboxylic acid is required, then the third fit shown in Table 1 would provide a more accurate value, but as a general predictive tool with a minimum of empirical parameters, the fourth fit (including all the dicarboxylic acids) should be more appropriate. 3.2. Phosphonic Acids. Our data set for the phosphonic acids (48 compounds, 90 pKa values) also included phosphoric acid and a few phosphates (see Table S2 in the Supporting Information). There are two pKa values among the phosphonic acids for which we found two different literature values which were fairly wide apart: pKa2 for methylphosphoric acid (7.2930 and 7.7631) and pKa2 for (3-chlorophenyl)phosphonic acid (6.6531 and 7.1032) (Table 3). On the basis of our linear regression fits, we have selected the first value shown in each case. The linear regression fit for the 90 P-OH pKa values gave R2 ) 0.98 with a mean absolute deviation of 0.28 pKa units. The fit is shown in Figure 2a and the corresponding linear equation is given in Table 1. The fit for the phosphonic acids is very similar to those for the monocarboxylic acids (see Table 1) so it is tempting to fit all the carboxylic and phosphonic acids (238 pKa values) together in a single fitting equation (Figure 2b). This gives a mean absolute deviation of 0.36 pKa units with R2 ) 0.94. Although this is a fairly good result (it offers a slight improvement over the mono- and dicarboxylic acids fit), in our opinion it is a significant deterioration from the phosphonic acid fit alone (the mean absolute deviation increases by 25% and the error range also increases) and so we cannot recommend it. 3.3. Phenols and Alcohols. Alcohols are commonly very weak acids unless there is a strong electron withdrawing group attached (for example fluorine; perfluoro-alcohols are quite strong acids). Phenols are usually considerably more acidic than alcohols. We used 67 phenols (rather loosely defined as anything involving an -OH group attached directly to an aromatic ring; see Table S3 in the Supporting Information) to develop the fitting equations. The resulting linear regression fit had a mean absolute deviation of 0.26 pKa units, with R2 ) 0.96 (Figure 3a). Although some authors subdivide the phenols into orthoand meta-/para-substituted subsets,33 which further improves the separate linear regression fits, this is inappropriate for our purposes, which is to produce a good predictiVe equation for a wide range of compounds, generally using as few fitting equations as possible. (See the comments made when discussing the carboxylic acids). Consequently, we prefer to consider all phenols together as a single group. Note that there are at least five phenols in our training set for which we found two or more literature experimental pKa values that were fairly wide apart (see Table S3 in the Supporting Information): 3-(methylsulfonyl)phenol (8.7530 and 9.3331,34); 4-(methylsulfonyl)phenol (7.8331,34 and 8.2830); 2,4,5trichlorophenol (7.0711 and 6.7235); 2-hydroxypyrazine (7.2836 and 8.2537); and 8-hydroxyquinoline (9.71 and 10.3736). On the basis of our linear regression fits, we have selected the first value shown in each case. The differences are quite significant, up to almost 1 pKa unit, showing the reduced reliability of the experimental data compared to those of paper 1.20

J. Phys. Chem. A, Vol. 114, No. 1, 2010 435 For the alcohols (28 compounds) the linear fit has a correlation coefficient of 0.97 and a mean absolute deviation of 0.40 pKa units (Figure 3b). The phenol and alcohol fits are very similar (Table 1) and strongly suggests, as with the aliphatic and aromatic carboxylic acids, that the two can be combined. The resulting linear regression fit (95 pKa values, Table 4) for what we will term the alcoholic OH group has a mean absolute deviation of 0.30 pKa units with R2 ) 0.98. 3.4. Hydroxamic Acids and Oximes. Both of these compound types contain the -N-OH group. We have already noted that hydroxamic acids with the -NH-OH motif can lose a proton from either the nitrogen or the oxygen; the resulting anion can be stabilized in either case by proton transfer: -NH-O- T -N--OH and so it would not appear to matter too much as to where the initial proton was lost. Consequently we consider all the hydroxamic acids together as a single group. The 45 hydroxamic acids in our data set (Table S4 in the Supporting Information) can be further divided into those with the full -NH-OH motif (33) and N-substituted species (12; with either a methyl or a phenyl group) where loss of the -OH proton is the only possibility. If all 45 are considered together, the resulting fit is clearly anomalous: the slope is 0.0582 and R2 ) 0.08 (although the mean absolute deviation is not so bad at 0.50 pKa units; see Table 5). Fitting the N-substituted hydroxamic acids separately from the rest clearly improves matters: R2 for these (12) is 0.60 with a mean absolute deviation of 0.17 pKa units (Figure 4b), while for the remainder (33) R2 ) 0.61 with a mean absolute deviation of 0.30 pKa units (Figure 4 c). Despite the improvement, the hydroxamic acids still give the worst correlation coefficients (R2) of all the compounds considered here. The linear regression fit to the oximes (27 compounds, 29 pKa values) had a mean absolute deviation of 0.39 pKa units with R2 ) 0.95. The fitting equation (see Table 1) was sufficiently similar to that for the hydroxamic acids (the -NH-OH acids not the N-substituted compounds, despite the presence of the -OH proton only) to suggest that the two data sets be combined. The resulting fit had a mean absolute deviation of 0.33 pKa units with R2 ) 0.92 (Figure 4d). 3.5. Thiols. We obtained the linear regression fit for the thiols (-SH groups) from a set of 40 compounds (42 pKa values; see Table S3 in the Supporting Information). This had a mean absolute deviation of 0.36 pKa units with a correlation coefficient (R2) of 0.95 (Figure 5a). The resulting fit (see Table 1) is reasonably close to that for the phenols + alcohols (the alcoholic OH group) and combining both groups (alcoholic OH + SH) gives a very similar fit with a mean absolute deviation of 0.34 pKa units and R2 ) 0.98 (Figure 5b). Following our data fitting we propose five general groups among the organic acids that we have tested: (1) the -COOH group (carboxylic acids) (2) the -POOH group (phosphorus acids) (3) the alcoholic -OH group (alcohols and phenols); this group can also include the -SH group (thiols) (4) the -NHOH and dNOH group (hydroxamic acids and oximes) (5) the -NROH group (N-substituted hydroxamic acids) The mean absolute deviation is 0.40 pKa units or less for all five groups with a maximum error range of (1.5 pKa units, with the sole exception of cis-cyclohexane-1,2-dicarboxylic acid for which the predicted pKa1 value (2.65) is 1.69 pKa units too low. Any (untested) molecule containing one (or more) of these functional groups for which the pKa is required should be

436

J. Phys. Chem. A, Vol. 114, No. 1, 2010

Zhang et al.

TABLE 2: Experimental and Calculated pKa Values for Carboxylic Acids pKaexp

pKacal

acetic acid bromoacetic acid butoxyacetic acid 2-(butylthio)acetic acid chloroacetic acid cyanoacetic acid dichloroacetic acid ethoxyacetic acid fluoroacetic acid isopropoxyacetic acid 2-(isopropyithio)acetic acid methoxyacetic acid 2-(methylsulfonyl)acetic acid nitroacetic acid 2-oxoacetic acid (glyoxylic acid) phenylacetic acid phenoxyacetic acid 2-(phenylsulfonyl)acetic acid 2-(phenylthio)acetic acid propoxyacetic acid 2-thiocyanatoacetic acid tribromoacetic acid trichloroacetic acid trifluoroacetic acid pivalic acid

4.76 2.90 3.66 3.74 2.84 2.47 1.35 3.84 2.72 3.69 3.72 3.57 2.36 1.48 2.98 4.31 3.17 2.51 3.57 3.69 2.58 0.72 0.52 0.53 5.05

4.30 2.77 3.66 3.76 2.72 2.52 1.88 3.93 2.85 3.71 3.66 3.59 2.75 1.65 2.10 4.29 3.26 3.00 3.49 3.50 2.35 0.92 0.99 0.84 4.73

benzoic acid 3-allylbenzoic acid 3-bromobenzoic acid 4-bromobenzoic acid 4-butoxybenzoic acid 3-chlorobenzoic acid 4-chlorobenzoic acid 3-cyanobenzoic acid 4-cyanobenzoic acid 2-fluorobenzoic acid 3-fluorobenzoic acid 4-fluorobenzoic acid 2-hydroxybenzoic acid 3-hydroxybenzoic acid 4-hydroxybenzoic acid 2-methylbenzoic acid 3-methylbenzoic acid

4.19 4.32 3.81 3.97 4.53 3.83 3.98 3.60 3.55 3.27 3.86 4.15 2.98 4.08 4.58 3.98 4.25

4.03 4.05 3.63 3.80 4.39 3.66 3.83 3.41 3.40 3.41 3.71 3.92 3.10 3.98 4.32 4.48 4.11

fumaric acid fumaric acid glutaric acid glutaric acid maleic acid maleic acid malonic acid malonic acid ethylmanolic acid ethylmanolic acid methylmalonic acid methylmalonic acid dimethylmalonic acid dimethylmalonic acid methylbutenedioic acid citraconic acid ethylmethylmalonic acid ethylmethylmalonic acid n-propylmalonic acid n-propylmalonic acid oxalic acid oxalic acid phthalic acid phthalic acid succinic acid succinic acid trans-caronic acid trans-caronic acid cyclopentane-1,1-dicarboxylic acid cyclopentane-1,1-dicarboxylic acid cyclopentane-1,1-diacetic acid 1,1-cyclopentanediacetic acid cis-cyclohexane-1,2-dicarboxylic acid

3.10 4.44 4.42 5.41 1.83 6.07 2.87 5.69 2.99 5.83 3.07 5.87 3.17 6.06 2.55 6.15 2.86 6.41 2.97 5.84 1.25 4.27 3.05 4.89 4.21 5.63 3.82 5.32 3.23 6.08 3.90 6.77 4.44 6.77

3.32 4.40 4.46 4.88 2.19 5.42 1.72 6.34 2.64 6.31 3.45 5.50 3.67 5.54 1.89 6.39 3.88 5.44 2.41 6.59 1.75 4.59 1.77 6.37 4.22 4.92 4.24 5.28 2.60 6.80 4.87 5.33 2.76 7.16

compounds

∆pKa

compounds

Aliphatic Carboxylic Acids 0.46 trans-2-chlorobut-2-enecarboxylic acid 0.13 cis-2-methylcyclopropanecarboxylic acid 0.00 trans-2-methylcyclo-propanecarboxylic acid -0.02 butyric acid 0.12 formic acid -0.05 glycolic acid -0.53 mercaptoacetic acid -0.09 3-methyl-2-butenoic acid -0.13 cyclohexanecarboxylic acid -0.02 3-formylpropanoic acid 0.06 hexanoic acid (caproic acid) -0.02 pentanoic acid (valeric acid) -0.39 ketomethylvaleric acid -0.17 propenoic acid (acrylic acid) 0.88 propargylic acid 0.02 pyruvic acid -0.09 acetopyruvic acid -0.49 propanoic acid 0.08 2-bromopropanoic acid 0.19 3-bromo-propanoic acid 0.23 2-chloropropanoic acid -0.20 3-chloropropanoic acid -0.47 lactic acid -0.31 2-methylpropanoic acid 0.32 Aromatic Carboxylic Acids 0.16 4-methylbenzoic acid 0.27 2-nitrobenzoic acid 0.18 3-nitrobenzoic acid 0.17 4-nitrobenzoic acid 0.14 2,4-dinitrobenzoic acid 0.17 3-methoxybenzoic acid 0.15 4-methoxybenzoic acid 0.19 3-propoxybenzoic acid 0.15 4-propoxybenzoic acid -0.14 2,3,4,5-tetramethylbenzoic acid 0.15 2,4,6-trimethylbenzoic acid 0.23 2,4,6-trinitrobenzoic acid -0.12 protocatechuic acid 0.10 2-furancarboxylic acid 0.26 naphthalene-1-carboxylic acid -0.50 2-phosphonobenzoic acid 0.14 Dicarboxylic Acids -0.22 trans-cyclohexane-1,2-dicarboxylic acid 0.04 -0.04 cis-cyclohexane-1,3-dicarboxylic acid 0.53 -0.36 trans-cyclohexane-1,3-dicarboxylic acid 0.65 1.15 trans-4-cyclohexene-1,2-dicarboxylic acid -0.65 0.35 trans-cyclopentane-1,2-dicarboxylic acid -0.48 -0.38 cis-cyclopentane-1,3-dicarboxylic acid 0.37 -0.50 trans-cyclopentane-1,3-dicarboxylic acid 0.52 0.66 trans-cyclopropane-1,2-dicarboxylic acid -0.24 -1.02 threo-2-bromo-3-chloro-succinic acid 0.97 0.56 rac-2,3-dibromosuccinic acid -0.75 rac-2,3-dibromosuccinic acid -0.50 meso-2,3-dibromosuccinic acid -0.32 meso-2,3-dibromosuccinic acid 1.28 erythro-2-bromo-3-chlorosuccinic acid -1.48 -0.01 meso-2,3-dichlorosuccinic acid 0.71 meso-2,3-dichlorosuccinic acid -0.42 rac-2,3-dichlorosuccinic acid 0.04 rac-2,3-dichlorosuccinic acid 0.63 meso-2,3-dimethylsuccinic acid -0.72 meso-2,3-dimethylsuccinic acid -1.07 rac-2,3-dimethylsuccinic acid 1.44 rac-2,3-dimethylsuccinic acid 1.69 -0.39

pKaexp

pKacal

∆pKa

3.22 5.02 5.00 4.82 3.76 3.83 3.58 5.12 4.90 3.60 4.85 4.86 2.30 4.25 1.89 2.60 2.61 4.87 3.01 4.01 2.89 3.99 3.83 4.64

3.45 4.94 4.78 4.63 3.22 3.33 3.65 4.93 4.83 4.18 4.64 4.65 2.47 4.11 2.21 2.68 2.79 4.63 3.31 3.77 3.24 3.91 3.45 4.78

-0.23 0.08 0.22 0.19 0.54 0.50 -0.07 0.19 0.07 -0.58 0.21 0.21 -0.17 0.14 -0.32 -0.08 -0.18 0.24 -0.30 0.24 -0.35 0.08 0.38 -0.14

4.37 2.21 3.46 3.43 1.42 4.09 4.50 4.20 4.46 4.22 3.55 0.65 4.48 3.27 4.24 1.71

4.20 2.52 3.30 3.27 1.68 4.03 4.36 4.06 4.39 4.48 4.81 1.29 4.30 3.36 4.26 1.46

0.17 -0.31 0.16 0.16 -0.26 0.06 0.14 0.14 0.07 -0.26 -1.26 -0.64 0.18 -0.09 -0.02 0.25

4.30 5.93 4.10 5.46 4.31 5.73 3.95 5.81 4.14 5.85 4.26 5.51 4.32 5.42 3.80 5.13 1.46 2.77 1.42 3.24 1.51 2.71 1.38 3.09 1.49 2.97 1.43 2.81 3.77 5.94 3.94 6.20

4.81 5.05 4.48 4.94 4.55 5.25 4.21 5.09 2.63 6.71 4.45 4.91 4.77 4.64 3.92 4.85 2.44 3.13 2.46 3.13 2.46 3.15 2.42 3.13 2.29 2.88 2.38 3.09 4.72 5.37 4.28 5.54

-0.51 0.88 -0.38 0.52 -0.24 0.48 -0.26 0.72 1.51 -0.86 -0.19 0.60 -0.45 0.78 -0.12 0.28 -0.98 -0.36 -1.04 0.11 -0.95 -0.44 -1.04 -0.04 -0.80 0.09 -0.95 -0.28 -0.95 0.57 -0.34 0.66

Method for Predicting pKa Values

J. Phys. Chem. A, Vol. 114, No. 1, 2010 437

TABLE 3: Experimental and Calculated pKa Values for Phosphonic Acids compounds

pKa1exp

pKa1cal

∆pKa

pKa2exp

pKa2cal

∆pKa

phosphonic acid methylphosphonic acid (2-methylbut-2-yl)phosphonic acid isobutylphosphonic acid n-butylphosphonic acid tert-butylphosphonic acid ethylphosphonic acid (bromomethyl)phosphonic acid chloromethylphosphonic acid (dichloromethyl)phosphonic acid (hydroxymethyl)phosphonic acid n-propylphosphonic acid isopropylphosphonic acid phenylphosphonic acid (3-hydroxylphenyl)phosphonic acid (2-bromophenyl)phosphonic acid (2-chlorophenyl)phosphonic acid (2-fluorophenyl)phosphonic acid (2-methoxyphenyl)phosphonic acid (3-chloro-4-methoxyphenyl)phosphonic acid (2-methylphenyl)phosphonic acid (3-methylphenyl)phosphonic acid (4-methylphenyl)phosphonic acid (3-nitrophenyl)phosphonic acid (2-chloro-4-nitrophenyl)phosphonic acid (4-sulfamoylphenyl)phosphonic acid (2,6-dimethylphenyl)phosphonic acid 2-phosphonobenzoic acid (2-hydroxy-4-nitrophenyl)phosphonic acid (2-nitrophenyl)phosphonic acid (2-bromo-5-methylphenyl)phosphonic acid (2-methoxy-4-nitrophenyl)phosphonic acid (2,4-dimethylphenyl)phosphonic acid (3,4-dimethylphenyl)phosphonic acid (3-chlorophenyl)phosphonic acid (4-cyanophenyl)phosphonic acid (2,4-dinitrophenyl)phosphonic acid benzylphosphonic acid [difluoro(phenyl)methyl]phosphonic acid [fluoro(phenyl)methyl]phosphonic acid (phenylmethyl)phosphonothioic acid (4-chlorophenyl)phosphorate (2,4-dichlorophenyl)phosphorate (2,4,5-trichlorophenyl)phosphorate (2-phenylvinyl)phosphonic acid propylthiophosphonic acid phenyl phosphate phosphoric acid

1.50 2.48 2.88 2.70 2.59 2.79 2.45 1.14 1.51 1.14 1.91 2.49 2.66 1.83 1.78 1.64 1.63 1.64 2.16 2.25 2.10 1.88 1.84 1.30 1.12 1.42 2.39

1.74 2.56 2.98 2.79 2.79 2.81 2.77 1.24 1.43 0.30 1.94 2.82 2.89 2.15 2.09 1.73 1.73 1.60 1.82 1.95 2.24 2.22 2.16 1.15 1.05 1.56 2.72

-0.24 -0.08 -0.10 -0.09 -0.20 -0.02 -0.32 -0.10 0.08 0.84 -0.03 -0.33 -0.23 -0.32 -0.31 -0.09 -0.10 0.04 0.34 0.30 -0.14 -0.34 -0.32 0.15 0.07 -0.14 -0.33

1.22 1.45 1.81 1.53 2.17 2.04 1.55 1.27

0.63 1.52 1.80 1.59 2.53 2.44 1.93 1.72

0.59 -0.07 0.01 -0.06 -0.36 -0.40 -0.38 -0.45

2.30

2.38

-0.08

0.70 0.60 0.53 2.00

0.51 0.38 0.25 2.23

0.19 0.22 0.28 -0.23

1.46 1.90

0.67 1.45

0.79 0.45

6.79 7.29 8.96 8.43 8.19 8.88 8.05 6.52 6.59 5.61 7.15 8.18 8.44 7.07 7.03 7.00 6.98 6.80 7.77 6.70 7.68 7.44 7.33 6.27 6.14 6.38 8.62 3.78 5.39 6.74 7.15 6.96 8.07 7.76 6.65 6.79 6.38 7.55 5.70 6.60 6.41 5.89 5.76 5.47 7.10 6.64 6.29 6.71

6.48 7.36 8.95 8.25 7.93 8.62 7.92 6.63 7.64 6.43 7.01 7.86 8.29 7.19 7.13 7.06 7.04 6.95 7.56 7.27 7.71 7.31 7.30 6.79 6.55 6.61 7.70 4.32 4.98 7.17 7.07 6.95 7.65 7.27 6.78 6.48 6.57 7.90 6.39 6.08 6.32 5.51 5.40 5.11 6.92 6.35 5.67 6.33

0.31 -0.07 0.01 0.18 0.26 0.26 0.13 -0.11 -1.05 -0.82 0.14 0.32 0.15 -0.12 -0.10 -0.06 -0.06 -0.15 0.21 -0.57 -0.03 0.13 0.03 -0.52 -0.41 -0.23 0.92 -0.54 0.41 -0.43 0.08 0.01 0.42 0.49 -0.13 0.31 -0.19 -0.35 -0.69 0.52 0.09 0.38 0.36 0.36 0.18 0.29 0.62 0.38

appropriately classified, the energy of the neutral and the anion obtained using our computational protocol, and the pKa estimated using the appropriate linear equation (see Table 1). Figure 6 depicts the distribution of residuals (the difference between the calculated and experimental pKa values) for each of the five functional groups. As can be seen, the bulk of the 449 pKa values that we have computed (almost 70% of them) lie within ( 0.4 pKa units of the experimental reference value. 4. Comparison with Prior Work We have already20 compared our computational protocol with earlier work of Adam,38 Friesner and co-workers,10 and Klamt and co-workers.11 Using a single fitting equation, Klamt and co-workers achieve a mean absolute deviation of 0.37 pKa units with a maximum deviation of 1.26 for a set of 64 organic and inorganic acids. With goals similar to ours, Friesner and coworkers fitted a ∼200 molecule training set with a mean unsigned deviation of 0.41 pKa units. However, their data set,

containing both acids and bases, was divided into 19 different functional groups with several of the functional groups containing as few as four compounds, which we consider too few for a three-parameter fit (they have an extra fitting parameter compared to our situation: the ionic radius in the solvation model). Their maximum deviation of 3.0 pKa units was also twice as large as ours. (This occurred for tert-butyl alcohol, a molecule that was also in our database and for which our error was just 0.33 pKa units; however, the experimental value used by Friesner and co-workers is principally to blame for this large deviation.) Adam has achieved an impressive mean absolute error of only 0.105 pKa units with a maximum of 0.342 for 19 aliphatic carboxylic acids. However, aliphatic carboxylic acids are one of the easiest groups to fit, and Adam’s database is much smaller and less varied than ours. The computational protocol in all three of these papers10,11,38 is more expensive than the one we use. Adam uses DFT with the PW91 functional39 (which is similar in cost to OLYP) and

438

J. Phys. Chem. A, Vol. 114, No. 1, 2010

Zhang et al.

Figure 2. Phosphonic acids.

also uses COSMO with the 6-311+G** basis set, again as we do, but his geometries are fully optimized with this basis, a much more costly option than our 3-21G(d) optimization. Klamt and co-workers carry out a full DFT geometry optimization using the BP86 functional40 with a TZVP basis set, all within their proprietary COSMO-RS method (COSMO for real solvents);41 again this is formally more expensive than our approach, although the actual cost is reduced as they use the numerically less accurate RI-DFT approach.42 They also perform a Boltzmann conformational average over potentially several low-energy conformers if this is considered necessary. (The precise circumstances as to when this is done are not stated, but it clearly adds to the computational cost.) The Friesner group approach is probably the most expensive, although again the formal cost is reduced because of their methodology.43 They first optimize geometries in the gas phase at B3LYP/6-31G*, followed by single-point energies at the B3LYP/cc-pvtz(-f) level44 with the further addition to an already very large basis a set of diffuse functions at the reactive center (presumably the site from which the proton is removed). Solvent effects are incorporated using a self-consistent reaction field approach.45 They were originally hoping that this direct approach would yield pKa values with a mean absolute deviation of 0.5 pKa units or less, but it took a further three-parameter linear regression fit to achieve this goal. As in this work, all three of the above authors use a linear regression fit. A number of other groups use an absolute approach which, as we have already indicated,20 is computationally very demanding and, unless very high levels of theory are used, often gives worse results. Such methods cannot be routinely used for fast throughput screening, which is one of our aims. An example is a recent paper from Namazian and Halvani.46 Their method involves full geometry optimization in the gas phase at the B3LYP/6-31+G** level and expensive explicitly computed thermodynamic corrections, with solvation

Figure 3. Phenols + alcohols.

energies calculated using a polarizable continuum model (PCM).47 These authors obtain an average error of 0.5 pKa units for a set of 66 carboxylic acids in aqueous solution. This is an excellent result for a direct method but is significantly worse than our mean absolute deviation (Table 1) of 0.20 pKa units over 82 monocarboxylic acids at far less computational cost. More typical is the earlier work of Silva and co-workers.9 Their model involved a full geometry optimization, in both the gas and aqueous phases, at HF/6-31+G**, with explicit thermodynamic corrections calculated in the gas phase (with frequencies scaled by 0.9181) and a PCM used in solution. There was also a “relaxation cycle” that involved additional singlepoint energies at the solution geometry in the gas phase and at the gas phase geometry in solution. They applied this model to a collection of 15 alcohols, thiols, and carboxylic acids, obtaining a rather poor mean absolute deviation of 1.73 pKa units, with a maximum error of 4.52. Unfortunately, results like this are all too common with direct approaches. Schmidt and Busch and Knapp8 also used a direct (nonfitted) approach to compute pKa’s for what they call a heterogeneous collection of 26 carboxylic acids, phenols. imides, pyridines, and imidazoles. They optimized geometries at the B3LYP/6-

Method for Predicting pKa Values

J. Phys. Chem. A, Vol. 114, No. 1, 2010 439

TABLE 4: Experimental and Calculated pKa Values for Phenols, Alcohols, and Thiols pKaexp

pKacal

∆pKa

compounds

pKaexp

pKacal

∆pKa

1-butanol ethanol 2,2-dichloroethanol 2-methoxyethanol trichloroethanol trifluoroethanol trinitroethanol methanol phenylmethanol 1-propanol 2-methylpropan-1-ol 2-propanol 1,4-butanediol 5-hydroxytetrazole

16.1 15.93 12.89 14.2 12.24 12.37 2.37 15.5 15.4 16.1 16.1 15.7 15.1 10.26

16.24 16.12 13.62 15.34 12.10 12.36 2.99 15.77 15.37 16.19 17.01 16.31 16.00 9.91

-0.14 -0.19 -0.73 -1.14 0.14 0.01 -0.62 -0.27 0.03 -0.09 -0.91 -0.61 -0.90 0.35

Alcohols 2,2,2-trifluoro-1-phenylethanol 2,2,2-trifluoro-1-(4-tolyl)-ethanol 2,2,2-trifluoro-1-(4-methoxyphenyl)ethanol 2,2,2-trifluoro-1-(3-nitrophenol)ethanol 2,2,3,3-tetrafluoro-1-propanol 3-(4-hydroxy-3-methoxyphenyl)-2-propen-1-ol 1,1,1-trifluoropropan-2-ol prop-2-en-1-ol tert-butyl alcohol 2-butanol 2-phenoxyethanol 1,2-ethanediol (ethylene glycol) 1,2-propanediol 1,3-propanediol

11.9 12.04 12.24 11.23 12.74 9.54 12.37 15.5 16.9 17.6 15.1 15.4 14.9 15.1

12.14 12.34 12.41 11.22 12.51 9.35 12.53 15.47 16.55 17.31 14.88 14.31 14.25 13.85

-0.24 -0.30 -0.17 0.01 0.23 0.19 -0.16 0.03 0.35 0.29 0.22 1.09 0.65 1.25

phenol 2-aminophenol 3-aminophenol 4-aminophenol 2-bromophenol 3-bromophenol 4-bromophenol 2-tert-butylphenol 3-tert-butylphenol 4-tert-butylphenol 2-chlorophenol 3-chlorophenol 4-chlorophenol 3-cyanophenol 4-cyanophenol 2,3-dichlorophenol 2,4-dichlorophenol 2,5-dichlorophenol 2,6-dichlorophenol 3,4-dichlorophenol 3,5-dichlorophenol 2,5-dimethylphenol 2,6-dimethylphenol 3,4-dimethylphenol 3,5-dimethylphenol 2,4-dinitrophenol 2,5-dinitrophenol 2,6-dinitrophenol 2-fluorophenol 3-fluorophenol 4-fluorophenol 2-methylphenol 3-methylphenol 4-methylphenol

9.98 9.96 9.96 10.46 8.45 8.87 9.35 10.62 10.12 10.23 8.53 8.88 9.38 8.57 7.95 7.71 7.89 7.51 6.81 8.62 8.18 10.41 10.6 10.36 10.19 4.12 5.2 3.73 8.7 9.29 9.81 10.32 10.09 10.27

9.83 10.40 10.52 11.31 8.38 8.61 9.06 10.84 10.31 10.30 8.43 8.76 9.17 8.20 7.28 7.62 7.78 7.37 6.93 8.25 7.75 10.58 10.61 10.47 10.37 4.08 5.22 4.24 8.56 8.91 9.62 10.30 10.09 10.24

0.15 -0.44 -0.56 -0.85 0.07 0.26 0.29 -0.22 -0.19 -0.07 0.10 0.12 0.21 0.37 0.67 0.09 0.11 0.14 -0.12 0.37 0.43 -0.17 -0.01 -0.11 -0.18 0.04 -0.02 -0.51 0.14 0.38 0.19 0.02 0.00 0.03

Phenols 2-methoxyphenol 3-methoxyphenol 4-methoxyphenol 3-(methylsulfonyl)phenol 4-methylsulfonylphenol 2-nitrophenol 3-nitrophenol 4-nitrophenol 2,4,5-trichlorophenol 4-hydroxyisoquinoline 5-hydroxyisoquinoline 6-hydroxyquinoline 7-hydroxyquinoline 4-amino-8-hydroxyquinoline 8-hydroxy-2-methylquinoline 2,3,4-trichlorophenol 2,4,6-trimethylphenol 2,4,6-trinitrophenol 2,3,4,6-tetrachlorophenol pentachlorophenol 4-hydroxybenzothiazole 6-hydroxyisoquinoline 7-hydroxyisoquinoline 8-hydroxyisoquinoline 2-hydroxypyridine 2-hydroxypyrazine 3-hydroxyquinoline 5-hydroxyquinoline 8-hydroxyquinoline 5-amino-8-hydroxyquinoline 8-hydroxy-4-methylquinoline 4-chloro-3,5-dimetthylphenol dichloroxylenol

9.98 9.65 10.24 8.75 7.83 7.22 8.36 7.14 7.07 8.7 8.47 8.9 8.87 10.71 10.16 7.1 10.89 0.37 5.62 4.9 8.85 9.17 8.9 8.42 11.99 7.28 8.08 8.56 9.71 11.24 9.99 9.7 8.28

10.31 9.85 10.49 8.18 7.61 7.04 7.90 5.74 6.85 8.15 8.39 8.71 8.63 10.60 10.13 7.12 11.05 1.72 5.65 5.07 8.83 8.11 8.55 8.07 11.44 7.03 8.19 8.41 9.84 11.55 10.10 9.73 8.52

-0.33 -0.20 -0.25 0.57 0.22 0.18 0.46 1.40 0.22 0.55 0.08 0.19 0.24 0.11 0.03 -0.02 -0.16 -1.35 -0.03 -0.17 0.02 1.06 0.35 0.35 0.55 0.25 -0.11 0.15 -0.13 -0.31 -0.11 -0.03 -0.24

thiophenol 2-methylbenzenethiol 3-methylbenzenethiol 4-methylbenzenethiol butane-1-thiol 2-methylbutane-2-thiol ethanethiol 2-ethoxyethanethiol methylthioglycolate ethyl-2-mercaptoacetate 2-mercaptoethanol methanethiol phenylmethanethiol 2-mercapto-2-methyl-1-propanol 3-mercaptopropane-1,2-diol prop-2-ene-1-thiol 2-methylpropane-2-thiol 2-propanethiol 2,3-dimercapto-1-propanol 2,3-dimercapto-1-propanol thioglycolic acid

6.62 6.64 6.66 6.82 10.66 11.22 10.61 9.38 8.08 7.95 9.72 10.33 9.43 9.85 9.46 9.96 11 10.86 8.62 10.57 10.4

7.05 7.59 7.22 7.28 10.12 10.61 10.01 9.40 8.04 7.65 9.38 9.76 9.49 10.34 8.98 9.12 10.53 10.20 7.81 10.54 10.52

-0.43 -0.95 -0.56 -0.46 0.54 0.61 0.60 -0.02 0.04 0.30 0.34 0.57 -0.06 -0.49 0.48 0.84 0.47 0.66 0.81 0.03 -0.12

Thiols 1,2-ethanedithiol 1,2-ethanedithiol 2-pyridylmethanethiol prop-1-ene-2-thiol cyclohexanethiol 5-mercaptouracil 3-methoxybenzenethiol 4-methoxybenzenethiol 3-chlorobenzenethiol 4-chlorobenzenethiol 4-bromobenzenethiol 4-nitrobenzenethiol 4-acetylbenzenethiol 3-nitrobenzenethiol trifluoroethanethiol 3-thiopropionic acid 3-mercaptopropan-1-ol methyl 3-mercaptopropanoate butane-1,4-dithiol 1,3-bis(2-mercaptoethyl)urea 2,3-dimercaptopropan-2-ol

9.05 10.56 8.82 7.86 10.69 5.3 6.39 6.78 5.78 6.14 6.02 4.72 5.33 5.24 7.30 10.27 10.19 9.33 9.88 9.26 9.04

9.23 10.90 9.15 7.74 10.26 6.11 7.07 7.39 6.35 6.58 6.52 4.27 5.69 5.70 7.57 10.55 9.82 9.46 9.71 9.34 8.51

-0.18 -0.34 -0.33 0.12 0.43 -0.81 -0.68 -0.61 -0.57 -0.44 -0.50 0.45 -0.36 -0.46 -0.27 -0.28 0.37 -0.13 0.17 -0.08 0.53

compounds

31G** level and computed single-point energies (using both B3LYP and the Becke half and half48 density functional) with the massive cc-pvqz basis set. Thermodynamic corrections were included, computed (like the optimized geometries) at the

B3LYP/6-31G** level. Their solvation model requires the estimation of atomic charges derived from a fit to the electrostatic potential and the solution of the Poisson equation (see ref 8 for details). Their best results, using the BHalfandHalf

440

J. Phys. Chem. A, Vol. 114, No. 1, 2010

Zhang et al.

TABLE 5: Experimental and Calculated pKa Values for Hydroxamic Acids and Oximes compounds ethanedial dioxime 1-chloro-ethanedial dioxime ethanediamide dioxime 2-oxo-2-(3-pyridyl)ethanal oxime 2-oxo-2-(4-pyridyl)ethanal oxime 1,2-cyclopentanedione dioxime 3-pentanone oxime 1,2-cyclohexanedione dioxime 1,2-cyclohexanedione dioxime trifluoroacetaldehyde oxime hexafluoroacetone oxime 1,2-benzoquinone dioxime 1-nitroethanal oxime 2,3-butanedione dioxime 2,3-butanedione dioxime benzohydroxamic acid 2-aminobenzohydroxamic acid 4-aminobenzohydroxamic acid 4-chlorobenzohydroxamic acid 2-chlorobenzohydroxamic acid 4-cyanobenzohydroxamic acid 2-fluorobenzohydroxamic acid 4-methoxybenzohydroxamic acid 4-fluorobenzohydroxamic acid 3-nitrobenzohydroxamic acid 4-nitrobenzohydroxamic acid salicylohydroxamic acid acetohydroxamic acid phenylacetohydroxamic acid n-butyrohydroxamic acid 2-pyrazinecarbohydroxamic acid pyrimide-2-carbohydroxamic acid N-phenyl-3-pyridinecarbohydroxamic acid N-phenyl acetohydroxamic acid N-methylbenzohydroxamic acid N-methyl-(4-cholorobenz)hydroxamic acid N-methyl-(4-methoxybenzo)hydroxamic acid N-methyl-(4-trifluoromethylbenzo)hydroxamic acid

pKaexp pKacal ∆pKa

compounds

8.88 8.35 11.37 7.80 7.80 10.10 12.60 10.68 12.16 8.90 6.00 6.93 7.40 10.42 12.05

8.94 7.72 11.45 7.77 7.57 10.20 11.87 10.26 12.57 8.26 6.48 7.54 7.10 10.25 12.55

Oximes -0.06 2-oxo-2-phenylethanal oxime 0.63 1-nitrosoethanal oxime -0.08 1,4-benzoquinone dioxime 1-O-methyl ester 0.03 2-bromo-1,4-benzoquinone dioxime 1-O-methyl ester 0.23 3-bromo-1,4-benzoquinone dioxime 1-O-methyl ester -0.10 2-chloro-1,4-benzoquinone dioxime 1-O-methyl ester 0.73 3-chloro-1,4-benzoquinone dioxime 1-O-methyl ester 0.42 1,4-benzoquinone-1-oxime -0.41 2-methyl-1,4-benzoquinone-1-oxime 0.64 3-methyl-1,4-benzoquinone-1-oxime -0.48 2-bromo-1,4-benzoquinone-1-oxime -0.61 3-bromo-1,4-benzoquinone-1-oxime 0.30 2-chloro-1,4-benzoquinone-1-oxime 0.17 3-chloro-1,4-benzoquinone-1-oxime -0.50

8.75 9.17 9.42 8.58 7.85 8.26 8.00 8.99 8.81 8.40 8.12 7.43 9.46 9.23 9.47 8.10 7.88 8.00 8.34 8.28 8.17 8.64 8.06

8.74 8.87 9.53 8.51 8.79 8.01 8.72 9.05 8.62 8.05 7.80 7.77 9.25 8.86 9.43 8.19 8.21 7.68 8.53 8.29 8.22 8.38 8.10

Hydroxamic Acids 0.01 formohydroxamic acid 0.30 hexanohydroxamic acid -0.11 1-naphthohydroxamic acid 0.07 nicotinhydroxamic acid -0.94 isonicotinhydroxamic acid 0.25 2-hydroxypropanehydroxamic acid -0.72 3-hydroxy-2-phenylpropanehydroxamic acid -0.06 pentanohydroxamic acid 0.19 picolinoylhydroxamic acid 0.35 propionohydroxamic acid 0.32 2-nitrobenzohydroxamic acid -0.34 2-methoxybenzohydroxamic acid 0.21 4-methylbenzohydroxamic acid 0.37 2-methylbenzohydroxamic acid 0.04 4-bromobenzohydroxamic acid -0.09 amino aceto (glycine)hydroxamic acid -0.33 0.32 N-methyl-(4-nitrobenzo) hydroxamic acid -0.19 N-methyl-(4-methylbenz)hydroxamic acid -0.01 N-methyl-(3-methylbenz)hydroxamic acid -0.05 N-methyl-(3-methoxy-4-methylbenz)hydroxamic acid 0.26 N-methyl-(3-nitro-4-chlorobenz)hydroxamic acid -0.04 N-methyl-(3,5-dinitrobenz)hydroxamic acid

density functional for the final single-point energies, had a good mean absolute deviation (0.44 pKa units) and a maximum deviation of 0.96, but the method is very expensive and their data set is rather small. Finally,wementionsomeQSPR(quantitativestructure-property relationship) work. This is a popular method that attempts to find a direct correlation between so-called molecular descriptors, often based entirely on structure, and the physical property of interest, in this case the pKa. Once the relationship has been determined, applying it to a new compound usually requires very little calculation and the quality of the results obtained, i.e., the error in the predicted pKa value, is often as good as or better than with far more computationally demanding techniques. In fact, we are of the opinion that as a compromise between accuracy and effort, a good QSPR method is difficult to beat. A problem with QSPR and similar methods is that they usually employ a large number of descriptors and fitted parameters, and thus they work partly by simply memorizing molecular fragments. This makes predictions for novel structural motifs that are not in the fitted database somewhat uncertain. A combination of first-principles based methods with QSPR-like descriptors appears ideal. Tehan and co-workers33,49 extracted sets of compounds from a large database, optimized their geometries at the semiempirical AM1 level, and attempted to find a correlation between various properties derived from frontier molecular orbital theory50 and the experimental pKa. These properties (termed “quantum mechanical descriptors”), nearly all functions of the molecular

pKaexp pKacal ∆pKa 8.25 5.81 9.20 8.70 8.60 8.70 8.60 6.35 6.95 6.90 5.70 5.55 5.70 5.60

8.07 6.14 8.51 7.89 7.86 7.95 7.95 6.55 6.68 6.91 5.87 5.97 5.97 6.04

0.18 -0.33 0.69 0.81 0.74 0.75 0.65 -0.20 0.27 -0.01 -0.17 -0.42 -0.27 -0.44

8.65 9.76 7.75 8.30 7.94 9.35 9.00 9.37 8.85 9.52 8.20 8.90 9.05 8.55 8.68 9.32

8.75 9.44 8.94 8.26 8.06 9.45 9.10 9.43 9.31 9.53 9.00 9.24 8.83 9.02 8.45 9.23

-0.10 0.32 -1.19 0.04 -0.12 -0.10 -0.10 -0.06 -0.46 -0.01 -0.80 -0.34 0.22 -0.47 0.23 0.09

7.91 8.46 8.52 8.41 7.78 7.41

7.76 8.34 8.33 8.36 7.98 7.88

0.15 0.12 0.19 0.05 -0.20 -0.47

orbital coefficients, included the electrophilic frontier electron density, the nucleophilic frontier electron density, the electrophilic superdelocalizability, the nucleophilic superdelocalizability, the radical superdelocalizability, and the atom selfpolarizability. Linear regression fits were attempted between one or more of these descriptors and the experimental pKa. Fitted data sets included 175 phenols, 99 benzoic acids, 185 aliphatic carboxylic acids, 55 anilines, and 77 amines. Fitting each data set “as is” did not always provide particularly good correlations, but further subdividing, for example, splitting the phenols into meta/para-substituted, ortho-hydrogen-bonded and the remaining ortho-substituted compounds, usually produced good fits, often using only one descriptor. The mean absolute deviation for a single-descriptor (the electrophilic superdelocalizability) fit to the 58 meta/para substituted phenols was 0.23 pKa units with a maximum deviation of under one unit. A two-descriptor fit to 143 aliphatic carboxylic acids (excluding 42 amino acids) resulted in a mean absolute deviation of 0.33 pKa units (cf., a mean absolute deviation of 0.39 pKa units over 149 carboxylic acids in this work, but this includes aliphatic, aromatic, and dicarboxylic acids (see Table 1, eq 4). However, the maximum error in refs 33 and 49 was over 3 pKa units and the correlation coefficient (average value of R2) was only 0.67. Crippen and co-workers51 used a large number of structural motifs (which they term SMARTS strings) such as any aromatic atom with a halogen substituent, any aromatic atom connected to a non-π-bonded carbon atom, a noncarbon aliphatic atom meta to a hydroxyl group and many, many more to try to find

Method for Predicting pKa Values

J. Phys. Chem. A, Vol. 114, No. 1, 2010 441

Figure 4. Hydroxamic acids + oximes.

Figure 6. Distribution of residuals for the calculated organic acids.

Figure 5. Thiols.

correlations among a set of 1881 unique monoprotic compounds involving 1088 ionizable oxygen, 33 ionizable sulfur, and 760 ionizable nitrogen atoms, respectively. After a fairly complicated (and expensive) refinement procedure, their final “fit” involved 169 SMARTS strings and could accurately predit the pKa in their data set with an R2 of 0.94 and an rms error of 0.68.

Lee and co-workers52 use a program called SPARC (Sparc Performs Automated Reasoning in Chemistry) to predict pKa values for a set of 123 compounds, including many known drugs. They obtained an R2 of 0.92 and an rms error of 0.78 pKa units. For a larger set of 537 compounds from an internal data set at Pfizer Global Research the same approach gave R2 ) 0.80 and an rms error of 1.05 pKa units. The SPARC computational approach is based on blending well-known established methods such as structure-activity relations, linear free energy relations, and semiempirical molecular orbital theory. 5. Conclusions We have developed linear regression equations based on a specific theoretical protocol20 to predict the acid dissociation constants (pKa) for a number of different acidic functional groups: (1) the -COOH group (carboxylic acids); (2) the

442

J. Phys. Chem. A, Vol. 114, No. 1, 2010

-POOH group (phosphorus acids); (3) the alcoholic -OH group (alcohols and phenols; also includes the “alcoholic” -SH group (thiols)); (4) the -NOH group (hydroxamic acids and oximes); and (5) the -NROH group (N-substituted hydroxamic acids). Within each group the mean absolute (unsigned) deviation between the theoretically calculated acidity constants using our linear regression equations and the experimental reference value is 0.4 pKa units or less, with a maximum error range of (1.5 pKa units. There is a single exception that falls outside this error range: cis-cyclohexane-1,2-dicarboxylic acid, which has an error of -1.69 pKa units. Over 68% of our computed pKa values lie within (0.4 pKa units of the experimental reference value, and over 91% of them lie within (0.8 pKa units. Our results compare very favorably with similar models developed by others which have far greater computational cost. Acknowledgment. This work has been supported by the National Science Foundation (Grant No. CHE-0515922), and the Mildred B. Cooper Chair at the University of Arkansas. Funds for the acquisition of the Red Diamond supercomputer were provided in part by the National Science Foundation under grant No. MRI-0421099. Supporting Information Available: Tables S1-S4 of the Supporting Information contain the names, structures, experimental pKa values (with references), and calculated energies of the protonated and deprotonated forms of these compounds, and deprotonation energies for all molecules discussed in the text. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Vo¨lgyi, G.; Ruiz, R.; Box, K.; Comer, J.; Bosch, E.; Taka´cs-Nova´k, K. Anal. Chim. Acta 2007, 583, 418. (2) Xie, Y.; Jiang, Y.; Ben-Amotz, D. Anal. Biochem. 2005, 343, 223. (3) Sjoeberg, H.; Karami, K.; Beronius, P.; Sundeloef, L.-O. Int. J. Pharm. 1996, 141, 63. (4) Popov, K.; Ronkkomaki, H.; Lajunen, L. H. J. J. Pure Appl. Chem. 2006, 78, 663. (5) Delatour, C.; Leclercq, L. Rapid Commun. Mass. Spectrom. 2005, 19, 1359. (6) Plasson, R.; Cottet, H. Anal. Chem. 2006, 78, 5394; 2007, 79, 3020. (7) Tajc, S. G.; Tolbert, B. S.; Basavappa, R.; Miller, B. L. J. Am. Chem. Soc. 2004, 126, 10508. (8) Schmidt am Busch, M; Knapp, E.-W. ChemPhysChem 2004, 5, 1513. (9) da Silva, C. O.; da Silva, E. C.; Nascimento, M. A. C. J. Phys. Chem. A 2000, 104, 2402. (10) Kliciæ, J. J.; Freisner, R. A.; Liu, S. Y.; Guida, W. C. J. Phys. Chem. A 2002, 106, 1327. (11) Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M. E. J. Phys. Chem. A 2003, 107, 9380. (12) Hollingsworth, C. A.; Seybold, P. G.; Hadad, C. M. Int. J. Quantum Chem. 2002, 90, 1396. (13) Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. J. Am. Chem. Soc. 2002, 124, 6421. (14) Jang, Y. H.; Sowers, L. C.; Cagin, T.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 274. (15) Murlowska, K.; Sadlej-Sosnowska, N. J. Phys. Chem. A 2005, 109, 5590. (16) Satchell, J. F.; Smith, B. J. Phys. Chem. Chem. Phys. 2002, 4, 4314. (17) Kaminski, G. A. J. Phys. Chem. B 2005, 109, 5884.

Zhang et al. (18) Han, J.; Tao, F-M J. Phys. Chem. A 2006, 110, 257. (19) Eckert, F.; Klamt, A. J. Comput. Chem. 2006, 27, 11. (20) Zhang, S.; Baker, J.; Pulay, P. J. Phys. Chem. A, DOI: 10.1021/ jp9067069. (21) Klamt, A.; Schu¨u¨rmann, G. J. Chem. Soc., Perkin Trans 2 1993, 799. (22) Hoe, W.-M.; Cohen, A. J.; Handy, N. C. Chem. Phys. Lett. 2001, 341, 319. (23) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650. (24) PQS, version 3.3; Parallel Quantum Solutions: 2013 Green Acres Road, Suite A, Fayetteville, AR 72703, U.S.A.; see http://www.pqschem.com (25) Red Diamond was supported by the National Science Foundation under grant number 0421099. (26) Chambers, R. W. Prog. Nucl. Acids Res. Mol. Biol. 1965, 5, 349. (27) Markman, A. L.; Ostrobrod, G. B.; Abdusalamova, E. K. IzV. Vyssh. Ucheb. ZaVed., Khim. Khim. Tekhnol. 1969, 12, 97. (28) PCModel v 8.5; Serena Software: Bloomington, IN, 2003. (29) Saunders, M. Houk, K. N.; Wu, Y-D; Still, W. C.; Lipton, M.; Chang, G.; Guida, W. J. Am. Chem. Soc. 1990, 112, 1419. (30) Serjeant, E. P.; Dempsey, B. Ionisation constants of organic acids in aqueous solution; IUPAC Chemistry Data Series No. 23; Pergamon Press: Oxford, U.K., 1979; see also references therein. (31) Kortum, G.; Vogel, W.; Andrussow, K. Dissociation constants of organic acids in aqueous solution; Butterworths: London, 1961; see also references therein. (32) Nagarajan, K.; Shelly, K. P.; Perkins, R. P.; Stewart, R. Can. J. Chem. 1987, 65, 1729. (33) Tehan, B, G.; Lloyd, E. J.; Wong, M. G.; Pitt, W. R.; Montana, J. G.; Manallack, D. T.; Gancia, E. Quant. Struct.-Act. Relat. 2002, 21, 457. (34) Brown, H. C.; McDaniel, H.; Hafliger, O. In Determination of organic structure by physical methods; Braude, E. A., Nachod, F. C., Eds.; Academic Press: New York, 1955; pp 567-662. (35) Ugland, K.; Lundanes, E.; Greibrokk, T. J. Chromatogr. 1981, 213, 83. (36) Perrin, D. D. Dissociation constants of organic bases in aqueous solution supplement; Butterworths: London, 1972; see also references therein. (37) Perrin, D. D. Dissociation constants of organic bases in aqueous solution; Butterworths: London, 1965; see also references therein. (38) Adam, K. R. J. Phys. Chem. A 2002, 106, 11963. (39) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. (40) Becke, A. D. Phys. ReV. A 1988, 38, 3098. Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (41) Klamt, A. COSMO-RS: From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design; Elsevier: Amsterdam, 2005. (42) Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 242, 652. (43) Friesner, R. A.; Murphy, R. B.; Beachy, M. D.; Ringnalda, M. N.; Pollard, W. T.; Dunietz, B. D.; Cao, Y. J. Phys. Chem. A 1999, 103, 1913. (44) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (45) Cortis, C. M.; Friesner, R. A. J. Comput. Chem. 1997, 18, 1570, 1591. (46) Namazian, M.; Halvani, S. J. Chem. Thermodyn. 2006, 38, 1495. (47) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Chem. Phys. Lett. 1996, 255, 327. (48) Becke, A. D. Phys. ReV. 1988, 38, 3098. (49) Tehan, B. G.; Lloyd, E. J.; Wong, M. G.; Pitt, W. R.; Gancia, E.; Manallack, D. T. Quant. Struct.-Act. Relat. 2002, 21, 473. (50) Fukui, K.; Yonezawa, T.; Nagata, C. Bull. Chem. Soc. Jpn. 1954, 27, 423. (51) Lee, A. C.; Yu, J.; Crippen, G. M. J. Chem. Inf. Model 2008, 48, 2042. (52) Lee, P. H.; Ayyampalayam, S. N.; Carreira, L. A.; Shalaeva, M.; Bhattachar, S.; Coselmon, R.; Poole, S.; Gifford, E.; Lombardo, F. Mol. Pharmacol. 2007, 4, 498.

JP9067087