Letter pubs.acs.org/JPCL
Resolution of a Challenge for Solvation Modeling: Calculation of Dicarboxylic Acid Dissociation Constants Using Mixed Discrete− Continuum Solvation Models Aleksandr V. Marenich, Wendu Ding, Christopher J. Cramer,* and Donald G. Truhlar* Department of Chemistry and Supercomputing Institute, University of Minnesota, 207 Pleasant Street Southeast, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *
ABSTRACT: First and second dissociation constants (pKa values) of oxalic acid, malonic acid, and adipic acid were computed by using a number of theoretical protocols based on density functional theory and using both continuum solvation models and mixed discrete−continuum solvation models. We show that fully implicit solvation models (in which the entire solvent is represented by a dielectric continuum) fail badly for dicarboxylic acids with mean unsigned errors (averaged over six pKa values) of 2.4−9.0 log units, depending on the particular implicit model used. The use of water−solute clusters and accounting for multiple conformations in solution significantly improve the performance of both generalized Born solvation models and models that solve the nonhomogeneous dielectric Poisson equation for bulk electrostatics. The four most successful models have mean unsigned errors of only 0.6−0.8 log units. SECTION: Molecular Structure, Quantum Chemistry, and General Theory An accurate prediction of the pKa values of polyprotic acids in aqueous solution is a long-standing and difficult challenge.1−3 We shall present calculations showing a dramatic failure of all available continuum solvation models for this important problem. Then, we shall show how explicit ensemble averaging over explicitly hydrated solutes surrounded by a continuum corrects the failure of models based solely on a continuum. The calculations are presented for three cases, oxalic acid (C2H2O4), malonic acid (C3H4O4), and adipic acid (C6H10O4) in water. Theory. The pKa value of a generic acid HA that dissociates into its conjugate base A− and a proton H+ in aqueous solution (aq) is determined by the following equation ΔGa°(aq) pK a = 2.3026RT
G°(M, aq) = G°(M·nH 2O, aq) − nG°(H 2O, aq) − nRT ln(55.34)
(3)
where the notation M·nH2O denotes a supersolute with n water molecules (nwat) added explicitly to the molecule M. The last term in eq 3 is a concentration term equal to 2.38 kcal/mol for n = 1 at 298 K (see details in refs 4 and 5). Scheme 1 illustrates the thermochemical cycle6 that relates ΔGa°(aq) to its gas-phase (g) counterpart, ΔG°a (g). One can use eq 3 to rewrite eq 2 as ΔGa°(aq) = G°(H+, aq) + G°(A−·nH 2O, aq) − G°(HA·nH 2O, aq)
(4)
(1)
where ΔG°a (aq) is the free energy of dissociation of HA defined as +
G°(X, aq) = −RT ln(
−
ΔGa°(aq) = G°(H , aq) + G°(A , aq) − G°(HA, aq)
∑ i ∈ {C}
(2)
e−Gi°(X,aq)/ RT ) (5)
where G°i(X,aq) is the free energy of conformer i. If there is only one conformer, then G°(X,aq) reduces to G°1 (X,aq). The free energy of conformer i is defined as
where G°(M,aq) is the standard-state free energy of species M in aqueous solution (M = H+, A−, or HA). Throughout the Letter, a small circle superscript refers to the 1 atm ideal-gas standard state for gases and the 1 M ideal-solution standard state for solutes. For a liquid-phase standard-state concentration of 1 M, the quantities G°(M,aq) obey the relation4 © 2012 American Chemical Society
−
The quantities G°(X,aq) (X = H , A ·nH2O, or HA·nH2O, where n ≥ 0) can be expressed as the thermodynamic average over a set {C} of conformers i according to +
Received: April 4, 2012 Accepted: May 6, 2012 Published: May 6, 2012 1437
dx.doi.org/10.1021/jz300416r | J. Phys. Chem. Lett. 2012, 3, 1437−1442
The Journal of Physical Chemistry Letters
Letter
Scheme 1. Thermochemical Cyclea
The energies Ei(X,g) used in eq 9 are calculated at the M062X/MG3S//M06-2X/MG3S level of theory. For the case of nwat = 0, we have also tested the CCSD(T)18/aug-ccpVTZ19,20//M06-2X/MG3S approach. The CCSD(T) singlepoint energy calculations (with frozen 1s orbitals on C and O) were carried out using the MOLPRO program.21 The quantity G°(X,aq) for an individual species X is calculated by eq 5 over the set of all geometrically unique and distinguishable conformers having Gi°(X,aq) values within 5 kcal/mol of the conformer found to have the lowest value of G°i (X,aq). To generate initial sets of trial conformers to be optimized by density functional theory (DFT), we carried out conformational searches with the MacroModel utility22,23 using the OPLS-200524,25 force field. The quantity ΔG*S,i in eq 7 was calculated for all species in eq 4 except H+ by using a dielectric continuum approximation with ε = 78.4 (the dielectric permittivity of liquid water at 298 K) at the M06-2X/MG3S gas-phase optimized geometry. For H+, we use the reference value of ΔG°S (H+) = −264.0 kcal/mol as reported by Tissandier et al.,26 which corresponds to ΔGS*(H+) = −265.9 kcal/mol. Continuum solvation calculations were carried out using M06-2X9,10 and B3LYP27−30 density functionals and the Hartree−Fock (HF) method and employed several basis sets (MG3S,11 cc-pVTZ,19 6-31+G(d,p),31 6-31G(d,p),31 6-31G(d),31 and SVP32). We tested the SMD model33 (solvation model based on density) by using Gaussian 09, the SM8 model34 (solvation model 8) and the SM8AD model35 (solvation model 8 with asymmetric descreening) by using a locally modified version of Gaussian 09, and the SMVLE36 (solvation model with volume and local electrostatics) as implemented in GAMESS37 (version of March 25, 2010, R2) as an extension of the SVPE38−40 (surface and volume polarization for electrostatics) model. The SVPE model itself was also tested, and we tested the conductorlike screening model (COSMO) protocol41,42 as implemented in Turbomole43 (version 6.3.1) and the Poisson−Boltzmann self-consistent reaction field solver44,45 (PBS) implemented in Jaguar.46 We also employed the pKa prediction module of Jaguar (PKA) utilizing a special multistep protocol that makes use of the PBS model and certain empirical corrections to the results of quantum mechanical calculations.46,47 The SMD, SM8, SM8AD, and SMVLE models (collectively referenced to as SMx models) were developed by the authors of this Letter previously,33−36 and these models were employed in the way described in the original publications. The other models tested here were employed with the default settings in the corresponding programs or by using the settings recommended in the programs’ manuals, unless noted otherwise. Examples of input files for solvation calculations are given in the Supporting Information (SI). The SMD model33 employs a self-consistent reaction field (SCRF) treatment of bulk electrostatics that involves solving the nonhomogeneous-dielectric Poisson equation (NPE) based on the continuous charge density of the solute. The Poisson equation is set up and solved by the IEF-PCM (integral equation formalism polarizable continuum model) computational algorithm.48 Nonbulk electrostatic solvation effects (which include both the nonelectrostatic effects and the deviation of the true electrostatics from the bulk electrostatics) are treated using the SMD cavity dispersion−solvent structure (CDS) protocols.33,49 The SMVLE36 and SVPE38−40 models can be interpreted as PCMs with exact treatments of volume polarization (which is
This relates the free energy of dissociation of HA in water, ΔG°a (aq), to the gas-phase free energy of dissociation, ΔGa°(g), using the standard-state solvation free energies of individual components, ΔG°S . a
Gi°(X, aq) = Gi°(X, g) + ΔGS,° i(X)
(6)
where the first term on the right side is the free energy of conformer i in the gas phase and the second term is the standard-state free energy of solvation,7 that is, of transfer from the gas phase at a solute partial pressure of 1 atm to a 1 M ° (X) is defined as aqueous solution at 298 K. The quantity ΔGS,i ΔGS,° i(X) = ΔGS,*i(X) + RT ln(RT /pV *)
(7)
where ΔGS,i * (X) is the fixed-concentration free energy of solvation equal to the free energy of transfer of the solute from the gas phase with a solute concentration of 1 mol/L to a 1 M solution. Using V* = 10−3 m3/mol (i.e., 1 L/mol) and p = 101325 Pa (i.e., 1 atm), we define the second term of eq 7 as the free-energy change of 1 mol of an ideal gas upon compression from 1 atm of pressure to 1 M concentration. At T = 298 K, this term equals ∼1.89 kcal/mol. The gas-phase free energy in eq 6 is defined as ° i(X, g) Gi°(X, g) = H0,° i(X, g) + Gtherm,
(8)
where the first term on the right side is the enthalpy of conformer i at 0 K defined by H0,° i(X, g) = Ei(X, g) + εZPE, i(X, g)
(9)
where the first term on the right side of eq 9 is the electronic energy (including nuclear repulsion) and the second term is the zero-point vibrational energy (ZPE). The second term of eq 8 is the thermal contribution to Gibbs free energy for the ith conformer in the gas phase. Computational Methods. The electronic energy and vibrational frequencies are calculated by Gaussian 098 using the M06-2X9,10/MG3S11 method at the corresponding gas-phase optimized geometry. Note that the MG3S basis set is identical to 6-311+G(3d2f,2df,2p) for elements with atomic numbers of 14 or lower, and it contains diffuse s and p subshells for all elements except hydrogen. The ZPE and gas-phase thermal free energy of all species except H+ are calculated from vibrational frequencies (scaled by a factor of 0.97012 to account for anharmonicity and for systematic errors in the electronic structure method) using the rigid rotor−harmonic oscillator approximation formulas,13 except that (after scaling) vibrational frequencies lower than 100 cm−1 are raised to 100 cm−1 as a way to correct for the well-known breakdown of the harmonic oscillator model for the free energies of low-frequency vibrational modes.14−16 In the case of H+ in the gas phase, we have H0,1 ° = 0 by definition, and Gtherm,1 ° is equal to −6.27 kcal/mol (see ref 17 and references therein). Because we use the harmonic oscillator formulas but scale the frequencies and raise the low frequencies to include some of the effects of anharmonicity, this may be called a quasiharmonic treatment. 1438
dx.doi.org/10.1021/jz300416r | J. Phys. Chem. Lett. 2012, 3, 1437−1442
The Journal of Physical Chemistry Letters
Letter
Table 1. pKa Values of Oxalic, Malonic, and Adipic Acidsa oxalic acid nwat
malonic acid nwat
pKa(1)
pKa(2)
1.25
3.81
0 4
−0.12 1.68
7.84 6.00
0 3
0
0.47
7.50
0
0 4
−1.27 0.66
3.18 4.29
0 3
0 4
−2.02 −0.23
0.27 3.85
0 3
0 4
−2.06 −0.39
−0.15 3.75
0 3
0 4
−4.79 −1.88
−4.26 −0.35
0 3
0 4
−3.85 −1.67
−2.69 1.64
0 3
0 4
−3.96 −1.38
−5.69 −0.32
0 3
0
−0.97
−18.65
0
0 4
−1.31 0.51
6.21 6.33
0 3
0 4
−3.72 −1.33
−1.37 2.16
0 3
0 4
−3.57 1.85
−4.21 −0.57
0 3
0
2.4
4.1
0
0
0.88
7.95
0
0
4.23
17.45
0
pKa(1)
adipic acid pKa(2)
experiment 2.85 5.70 SMD/M06-2X/MG3S 0.91 12.85 3.20 8.28 SMD/M06-2X/MG3Sb 2.12 12.27 SMD/M06-2X/cc-pVTZ −0.06 9.60 2.42 6.51 SMD/M06-2X/6-31G(d) −0.93 7.41 1.72 5.76 SMD/B3LYP/6-31G(d) −0.87 7.41 1.73 5.67 SM8/M06-2X/6-31+G(d,p) −3.00 3.86 0.02 2.28 SM8/M06-2X/6-31G(d) −1.72 4.67 0.96 3.56 SM8AD/M06-2X/6-31G(d) −1.56 1.50 1.82 1.42 SMVLE/HF/6-31+G(d) −1.69 −9.67 COSMO/B3LYP/SVP −0.72 13.70 1.91 8.58 COSMO/B3LYP/SVP, with SMD radii −2.98 6.27 −0.20 4.35 PBS/B3LYP/6-31G(d,p) −0.67 1.84 2.20 3.51 PKA 2.5 6.4 SVPE/B3LYP/6-31G(d) 1.39 14.89 SVPE/HF/6-31+G(d) 3.95 21.87
nwat
pKa(1)
pKa(2)
MUE
4.41
5.41
0 2
6.04 5.65
11.15 8.94
3.6 1.7
0
7.35
10.77
3.3
0 2
5.53 5.00
8.97 7.56
2.4 0.8
0 2
5.10 4.42
7.35 6.57
2.5 0.6
0 2
5.39 4.58
7.39 6.82
2.6 0.7
0 2
1.62 1.61
8.77 4.39
4.7 2.9
0 2
4.24 3.61
7.17 5.94
3.2 1.7
0 2
3.84 3.53
5.47 3.68
4.0 2.4
0
3.18
−2.92
9.0
0 2
5.51 5.16
12.21 9.96
4.1 2.1
0 2
3.42 3.19
5.74 4.96
3.0 1.7
0 2
3.42 6.03
2.75 1.26
4.0 2.3
0
5.7
5.7
0.7
0
7.73
12.63
4.3
0
9.25
18.37
8.6
a
The quantity nwat is the number of explicit water molecules added to the solute molecule. MUE is the mean unsigned error calculated over six pKa values relative to experimental reference. bUsing the Ei(X,g) values (eq 9) calculated by CCSD(T)/aug-cc-pVTZ//M06-2X/MG3S.
according to refs 44−46; as discussed elsewhere,50 this is very much in the spirit of the SMx models. The SM834 and SM8AD35 models are self-consistent reaction field continuum solvation models based on the generalized Born (GB) approximation51,52 for the bulk electrostatic contribution to the free energy of solvation, thereby representing the solute’s charge density in terms of point charges (a distributed monopole approximation), located at the nuclear positions. The SM8 model is defined using the effective Born radius formula proposed originally by Still et al.53 (with one parameter reoptimized54) using an unshielded Coulomb field for the electric displacement field induced by a point charge. The SM8AD model uses the asymmetric descreening
an additional polarization effect due to the solute’s density penetrating outside of the solute’s molecular cavity). The SMVLE model improves on the SVPE model by using a special term to account for the local electrostatic effect derived from the outward-directed normal electric field on the cavity surface and by using the CDS protocol to account for nonelectrostatic effects.36 The COSMO model41,42 in Turbomole43 solves the NPE in terms of screening surface charges obeying a conductor-like equation followed by scaling by a factor of (ε − 1)/(ε + 1/2), but it does not include any nonbulk electrostatic terms. The PBS model of the Jaguar program solves the NPE by its own unique combination of computational protocols,44−46 and it treats nonbulk electrostatic solvation effects in water 1439
dx.doi.org/10.1021/jz300416r | J. Phys. Chem. Lett. 2012, 3, 1437−1442
The Journal of Physical Chemistry Letters
Letter
Table 2. An Example of pKa Values Calculated Without Account For Multiple Conformersa,b oxalic acid
malonic acid
nwat
pKa(1)
pKa(2)
nwat
0 4
−2.45 1.18
0.27 4.74
0 3
pKa(1)
adipic acid pKa(2)
SMD/M06-2X/6-31G(d) −1.13 7.51 1.40 5.82
nwat
pKa(1)
pKa(2)
0 2
3.88 2.92
7.38 9.20
MUE 2.6 1.3
Using only the gas-phase global minimum for each species corresponding to the conformer found to have the lowest value of G°i (X,g). bSee footnote a to Table 1. a
algorithm of Grycuk55 rather than the Coulomb field approximation to treat dielectric descreening effects. For bulk electrostatic calculations, the SMVLE and SVPE models employ the solute’s dielectric cavity determined selfconsistently based on an isocontour of the solute electron density (we chose 0.001 atomic unit as a generally recommended isodensity contour value), whereas the remaining models employ the solute−solvent electrostatic boundary based on a set of intrinsic atomic Coulomb radii whose values differ between different models in general. We have tested mixed discrete−continuum solvation models (except for PKA, SMVLE, and SVPE) by adding nwat explicit water molecules to a bare solute molecule representing the neutral, the anion, and the dianion of the corresponding dicarboxylic acid, the maximum value of nwat being equal to 4, 3, and 2 for oxalic acid, malonic acid, and adipic acid, respectively. Table 1 shows the pKa(1) and pKa(2) values computed using eq 5 and the protocol described above in comparison with available experimental56 values. We see that all of the calculations (except the PKA ones) with no explicit water molecules have large errors. This is a crisis for continuum solvation modeling. However, we also see that the crisis is averted by the present mixed discrete− continuum method; adding explicit water molecules yields results substantially more accurate than those obtained in the case of nwat =0 across all of the solvation protocols tested. For example, the mean unsigned error (MUE) in the pKa values computed using the SMD/M06-2X/6-31G(d) protocol changes from 2.5 log units (nwat = 0) to 0.6 (nwat > 0) (Table 1). The use of the CCSD(T)/aug-cc-pVTZ//M06-2X/MG3S method to calculate the values of Ei(X,g) used in eq 9 improves the DFT results only a little on average (Table 1, footnote b). Overall, the SMD model outperforms any other model in Table 1, including those that also solve the NPE (namely, COSMO, PBS, SMVLE, and SVPE), except for the PKA prediction module of Jaguar, which shows comparable accuracy to SMD. We note though that the PKA model makes use of empirical corrections specific to certain functional groups in a molecule (in this case, the COO− moiety) to improve the ab initio computed pKa values, whereas the present approach does not require such add-on corrections. The SMD performance slightly worsens when the MG3S basis (containing diffuse functions) is employed in solvation calculations; this is possibly due to an increase in the outlying charge error associated with the use of a more diffuse basis set within the PCM approach.57 The PBS model in our calculations incorrectly predicts that pKa(2) < pKa(1) for oxalic and adipic acids. The trends predicted by the COSMO model are qualitatively correct, and the model’s average accuracy can be slightly improved by replacing its default values58 for intrinsic Coulomb radii with those for SMD,33 as shown in Table 1. The performance of the SMVLE model,
which was parametrized only for the HF/6-31+G(d) level of theory,36 is generally poor, especially for pKa(2), because the SMVLE/HF/6-31+G(d) solvation free energy of the corresponding dianion is too negative. The similar SVPE model performs better with the 6-31G(d) basis set. The tested GB models (SM8 and SM8AD) perform poorly for the pKa values of oxalic acid, but they yield more accurate results for adipic acid because the corresponding COO− moieties are more separated in the latter case. This is in agreement with the analysis of pKa shifts in dicarboxylic acids undertaken by Jayaram et al.,1 which focused on the dielectric shielding of charges by the solvent and the self-polarization terms arising within the GB approach. In the current study, the SM8AD model,35 which utilizes a new formula for the effective Born radius based on the asymmetric descreening algorithm, shows no improvement over the older SM8 model that uses the Coulomb field approximation to treat dielectric descreening effects. The use of a basis set with diffuse functions, for example, 6-31+G(d,p) versus 6-31G(d), in solvation calculations that employ the GB models does not improve the resulting pKa values (see Table 1). Table 2 shows the pKa values computed without averaging the quantity G°(X,aq) used in eq 2 over multiple conformers in solution by eq 5. In this calculation, we instead used only the global minimum gas-phase optimized conformer. This approach is traditionally used in the literature (see, for example, ref 59 and references therein), but it may be inaccurate because it neglects the additional entropy contributions that would come from considering additional conformations, especially when such contributions become important (the case of nwat > 0). A conformational averaging protocol for implicit solvation calculations on the pKa values of carboxylic acids was previously employed by Liptak and Shields,60 though they did not use mixed discrete−continuum models. The pKa(1) value of oxalic acid reported in ref 60 is somewhat different from ours (the case of nwat = 0 in Table 1) because of using different reference values for the solvation free energy of the proton in water, electronic structure methods, and continuum solvation models. The importance of including an ensemble average over multiple structures has also been emphasized recently for gas-phase thermochemistry.61 We have found that pure dielectric continuum models fail for dicarboxylic acids, which involve the screening of highly concentrated charges in close proximity to one another, and we expect this observation to be more general and to extend beyond the current demonstration set of dicarboxylic acids. We have shown that explicit consideration of a few first-shell solvent molecules can ameliorate this deficiency. We have noted previously the potential importance of including firstshell solvent molecules in cases of one concentrated charge center in a molecule,62 and the current study confirms its logical extension to multiply charged centers, where the effects are dramatic. 1440
dx.doi.org/10.1021/jz300416r | J. Phys. Chem. Lett. 2012, 3, 1437−1442
The Journal of Physical Chemistry Letters
■
Letter
(10) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (11) Lynch, B. J.; Zhao, Y.; Truhlar, D. G. Effectiveness of Diffuse Basis Functions for Calculating Relative Energies by Density Functional Theory. J. Phys. Chem. A 2003, 107, 1384−1388. (12) Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G. Computational Thermochemistry: Scale Factor Databases and Scale Factors for Vibrational Frequencies Obtained from Electronic Model Chemistries. J. Chem. Theory Comput. 2010, 6, 2872−2887. (13) Ochterski, J. W. Thermochemistry in Gaussian; Gaussian, Inc., 2000; http://www.gaussian.com/g_whitepap/thermo/thermo.pdf (accessed January 30, 2012). (14) Pratt, L. M.; Truhlar, D. G.; Cramer, C. J.; Kass, S. R.; Thompson, J. D.; Xidos, J. D. Aggregation of Alkyllithiums in Tetrahydrofuran. J. Org. Chem. 2007, 72, 2962−2966. (15) Zhao, Y.; Truhlar, D. G. Computational Characterization and Modeling of Buckyball Tweezers: Density Functional Study of Concave-Convex π···π Interactions. Phys. Chem. Chem. Phys. 2008, 10, 2813−2818. (16) Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Use of Solution-Phase Vibrational Frequencies in Continuum Models for the Free Energy of Solvation. J. Phys. Chem. B 2011, 115, 14556− 14562. (17) Zhan, C.-G.; Dixon, D. A. Absolute Hydration Free Energy of the Proton from First-Principles Electronic Structure Calculations. J. Phys. Chem. A 2001, 105, 11534−11540. (18) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A Fifth-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479−483. (19) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (20) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. Electron Affinities of the First-Row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (21) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A. MOLPRO, version 2010.1, A Package of Ab Initio Programs; 2010; Cardiff University: Cardiff, U.K. See http://www.molpro.net. (22) MacroModel, version 9.6; Schrodinger, LLC: New York, 2008. (23) Maestro, version 8.5.111 and MMshare, version 1.7.110; Schrodinger, LLC: New York, 2008. (24) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (25) Banks, J. L.; Beard, H. S.; Cao, Y.; Cho, A. E.; Damm, W.; Farid, R.; Felts, A. K.; Halgren, T. A.; Mainz, D. T.; Maple, J. R.; Murphy, R.; Philipp, D. M.; Repasky, M. P.; Zhang, L. Y.; Berne, B. J.; Friesner, R. A.; Gallicchio, E.; Levy, R. M. Integrated Modeling Program, Applied Chemical Theory (IMPACT). J. Comput. Chem. 2005, 26, 1752−1780. (26) Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R., Jr. The Proton’s Absolute Aqueous Enthalpy and Gibbs Free Energy of Solvation from Cluster-Ion Solvation Data. J. Phys. Chem. A 1998, 102, 7787−7794. (27) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100.
ASSOCIATED CONTENT
S Supporting Information *
An extended set of the computed pKa values and fixedconcentration free energies of solvation for studied species, examples of input files (Part I), and Cartesian coordinates for tested molecular structures (Part II). This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected] (C.J.C);
[email protected] (D.G.T.). Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported in part by the U.S. Army Research Laboratory under Grant W911NF09-100377, by the National Science Foundation under Grants CHE09-56776 and CHE0952054, and by an NSF/Lando fellowship to W.D. A portion of the research (MOLPRO calculations) was performed using EMSL, a national scientific user facility sponsored by the Department of Energy’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory, under Project gc34900.
■
REFERENCES
(1) Jayaram, B.; Liu, Y.; Beveridge, D. L. A Modification of the Generalized Born Theory for Improved Estimates of Solvation Energies and pK Shifts. J. Chem. Phys. 1998, 109, 1465−1471. (2) Smiechowski, M. Theoretical Calculation of pKa’s of Phosphoric (V) Acid in the Polarisable Continuum and Cluster-Continuum Models. J. Mol. Struct. 2009, 924−926, 170−174. (3) Lee, T. B.; McKee, M. L. Dependence of pKa on Solute Cavity for Diprotic and Triprotic Acids. Phys. Chem. Chem. Phys. 2011, 13, 10258−10269. (4) Bryantsev, V. S.; Diallo, M. S.; Goddard, W. A., III. Calculation of Solvation Free Energies of Charged Solutes Using Mixed Cluster/ Continuum Models. J. Phys. Chem. B 2008, 112, 9709−9719. (5) Pliego, J. R. Reply to Comment on: ‘Thermodynamic Cycles and the Calculation of pKa’ [Chem. Phys. Lett. 367 (2003) 145]. Chem. Phys. Lett. 2003, 381, 246−247. (6) Alongi, K. S.; Shields, G. C. Theoretical Calculations of Acid Dissociation Constants: A Review Article. Annu. Rep. Comput. Chem. 2010, 6, 113−138. (7) Ben-Naim, A. Solvation Thermodynamics; Plenum: New York, 1987; p 4. (8) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (9) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. 1441
dx.doi.org/10.1021/jz300416r | J. Phys. Chem. Lett. 2012, 3, 1437−1442
The Journal of Physical Chemistry Letters
Letter
Functional Theory and Self-Consistent Reaction Field Methods. J. Phys. Chem. A 2002, 106, 1327−1335. (48) Tomasi, J.; Mennucci, B.; Cancès, E. The IEF Version of the PCM Solvation Method: An Overview of a New Method Addressed to Study Molecular Solutes at the QM Ab Initio Level. J. Mol. Struct.: THEOCHEM 1999, 464, 211−226. (49) Liotard, D. A.; Hawkins, G. D.; Lynch, G. C.; Cramer, C. J.; Truhlar, D. G. Improved Methods for Semiempirical Solvation Models. J. Comput. Chem. 1995, 16, 422−440. (50) Cramer, C. J.; Truhlar, D. G. Reply to Comment on “A Universal Approach to Solvation Modeling”. Acc. Chem. Res. 2009, 42, 493−497. (51) Born, M. Volumen und Hydratationswärme der Ionen. Z. Phys. 1920, 1, 45−48. (52) Klopman, G. Solvations: A Semi-Empirical Procedure for Including Solvation in Quantum Mechanical Calculations of Large Molecules. Chem. Phys. Lett. 1967, 1, 200−202. (53) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129. (54) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. SM6: A Density Functional Theory Continuum Solvation Model for Calculating Aqueous Solvation Free Energies of Neutrals, Ions, and Solute− Water Clusters. J. Chem. Theory Comput. 2005, 1, 1133−1152. (55) Grycuk, T. Deficiency of the Coulomb-Field Approximation in the Generalized Born Model: An Improved Formula for Born Radii Evaluation. J. Chem. Phys. 2003, 119, 4817−4826. (56) CRC Handbook of Chemistry and Physics, 91st ed. (2010−2011); Haynes, W. M., Ed.; CRC Press: Boca Raton, FL, 2010. (57) Baldridge, K.; Klamt, A. First Principles Implementation of Solvent Effects without Outlying Charge Error. J. Chem. Phys. 1997, 106, 6622−6633. (58) Klamt, A.; Eckert, F. COSMO-RS: A Novel and Efficient Method for the A Priori Prediction of Thermophysical Data of Liquids. Fluid Phase Equilib. 2000, 172, 43−72. (59) 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. (60) Liptak, M. D.; Shields, G. C. Accurate pKa Calculations for Carboxylic Acids Using Complete Basis Set and Gaussian-n Models Combined with CPCM Continuum Solvation Methods. J. Am. Chem. Soc. 2001, 123, 7314−7319. (61) Zheng, J.; Yu, T.; Papajak, E.; Alecu, I. M.; Mielke, S. L.; Truhlar, D. G. Practical Methods for Including Torsional Anharmonicity in Thermochemical Calculations on Complex Molecules: The Internal-Coordinate Multi-Structural Approximation. Phys. Chem. Chem. Phys. 2011, 13, 10885−10907. (62) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Adding Explicit Solvent Molecules to Continuum Solvent Calculations for the Calculation of Aqueous Acid Dissociation Constants. J. Phys. Chem. A 2006, 110, 2493−2499.
(28) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle− Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (29) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (30) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (31) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (32) Schäfer, A.; Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571− 2577. (33) 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. (34) 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. (35) 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. (36) 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. (37) 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.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (38) Chipman, D. M. Charge Penetration in Dielectric Models of Solvation. J. Chem. Phys. 1997, 106, 10194−10206. (39) Zhan, C.-G.; Bentley, J.; Chipman, D. M. Volume Polarization in Reaction Field Theory. J. Chem. Phys. 1998, 108, 177−192. (40) Chipman, D. M. New Formulation and Implementation for Volume Polarization in Dielectric Continuum Theory. J. Chem. Phys. 2006, 124, 224111/1−224111/10. (41) Klamt, A.; Schüürmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (42) Schäfer, A.; Klamt, A.; Sattel, D.; Lohrenz, J. C. W.; Eckert, F. COSMO Implementation in TURBOMOLE: Extension of an Efficient Quantum Chemical Code Towards Liquid Systems. Phys. Chem. Chem. Phys. 2000, 2, 2187−2193. (43) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (44) Tannor, D. J.; Marten, B.; Murphy, R.; Friesner, R. A.; Sitkoff, D.; Nicholls, A.; Ringnalda, M.; Goddard, W. A., III; Honig, B. Accurate First Principles Calculation of Molecular Charge Distributions and Solvation Energies from Ab Initio Quantum Mechanics and Continuum Dielectric Theory. J. Am. Chem. Soc. 1994, 116, 11875− 11882. (45) Marten, B.; Kim, K.; Cortis, C.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.; Honig, B. New Model for Calculation of Solvation Free Energies: Correction of Self-Consistent Reaction Field Continuum Dielectric Theory for Short-Range Hydrogen-Bonding Effects. J. Phys. Chem. 1996, 100, 11775−11788. (46) Jaguar 7.7, release 107; Schrodinger, LLC: New York, NY, 2010. (47) Klicić, J. J.; Friesner, R. A.; Liu, S.-Y.; Guida, W. C. Accurate Prediction of Acidity Constants in Aqueous Solution via Density 1442
dx.doi.org/10.1021/jz300416r | J. Phys. Chem. Lett. 2012, 3, 1437−1442
