Modeling Flexible Molecules in Solution: A pKa Case Study - The

Jun 30, 2017 - Here we explore these issues using the calculation of pKas for a diverse set of amine-containing species as a case study. We show that ...
0 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCA

Modeling Flexible Molecules in Solution: A pKa Case Study Naomi L. Haworth,* Qinrui Wang, and Michelle L. Coote* ARC Centre of Excellence for Electromaterials Science, Research School of Chemistry, Australian National University, Canberra, ACT 2601, Australia S Supporting Information *

ABSTRACT: Continuum solvation models have been incredibly successful for the computationally efficient study of chemical reactions in solution. However, their development and application has generally been on focused on investigations of small, rigid molecules. Additional factors must be considered when studying large, flexible and multiply ionizable species. These include whether the use of thermocycle or entirely solution-phase approaches are more appropriate for the calculation of solution-phase free energies, which metrics can be used to reliably identify the conformation(s) adopted by flexible molecules in solution, and how errors due to inaccuracies in the prediction of low energy vibrational frequencies can be avoided. Here we explore these issues using the calculation of pKas for a diverse set of aminecontaining species as a case study. We show that thermocycle-based approaches should only be applied where there are relatively small structural changes between the gas- and solution-phase molecular geometries, and that these methods are generally not appropriate for conformational searching. Using gas- or solution-phase energies or gas-phase free energies can also lead to errors in the identification of the most stable molecular conformation(s). Scaling of low energy vibrational modes (i.e., use of the quasiharmonic oscillator approximation) is helpful, however care must be taken to ensure modes that change as part of the reaction are not disregarded. Entirely solution-phase approaches to the Gibbs free energy and hence pKa calculations were found to yield accurate pKa values for the amine test set studied when each charged site is complexed with an explicit water molecule and a proton exchange scheme is applied with an appropriately chosen reference acid.

■

INTRODUCTION The development of accurate, low-cost methods for modeling chemical behavior in solution is arguably one of the biggest ongoing challenges for computational chemistry. Although it would be desirable to describe solution-phase phenomena using an ensemble of explicit solvent molecules, for many purposes this is not computationally feasible. Instead, implicit solvation methods have been developed where the solvent is treated as a polarizable dielectric continuum.1−10 While these methods have proven to be highly successful for modeling solvent effects,11 it is important to recognize the development of approaches for their application to chemical problems have in most cases focused on test sets of small, rigid molecules. Many systems of chemical and pharmacological interest are neither small nor rigid, however. Here we explore important factors that must additionally be considered when modeling solution-phase properties of larger, more flexible species. Chief among these is the question of how to accurately explore the conformational space of flexible molecules. Different molecular conformations can have different stabilities and properties. For conformationally flexible species it is therefore necessary to identify the conformer(s) adopted by the majority of molecules in the relevant environment; for solution-phase chemistry, this means the conformation(s) with the lowest free energy in the relevant solvent. As solution-phase free energy © 2017 American Chemical Society

calculations can be computationally expensive, conformational searching is often performed in the gas phase12−14 or, if solutionphase calculations are performed, only energies are compared.12,13,15−17 However, significant conformational and structural changes (e.g., tautomerization) can occur between the gas and solution phases, particularly for larger and more complex species, and free energy trends do not always match those of the absolute energies. For example, intramolecular interactions such as hydrogen bonds and dispersion forces are energetically favored, particularly in the gas phase. However, this is often offset by a reduction in the associated entropy, and a corresponding increase in the free energy, in solution. Moreover, depending on the nature of the solvent and the hydrophobicity of solute, solvation also tends to favor more open conformations.18 A second important issue lies in how implicit solvation models are applied to chemical problems. Originally, these methods were designed to predict the Gibbs free energies of solvation of the relevant species. These solvation free energies were intended to be applied as corrections to accurate gas-phase free energies via a Hess’ law thermochemical cycle (thermocycle) to convert these into solution-phase values (see Scheme 1). This works well for Received: May 2, 2017 Revised: June 16, 2017 Published: June 30, 2017 5217

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A Scheme 1. Generic Thermocycle for the Calculation of the Solution-Phase Gibbs Free Energy of a Chemical Reaction Using the Gas-Phase Reaction Energy along with the Gibbs Free Energies of Solvation and the Associated Structural Relaxation of Each Species

Chart 1. Amine Test Set Used in this Work

small, rigid molecules, where there is little difference between the gas- and solution-phase geometries, and hence ΔGrelax (which cannot actually be calculated) can be approximated by ΔErelax. Difficulties may arise for larger systems, however, where significant geometric relaxation, accompanied by changes in vibrational frequencies and hence failure of ΔGrelax ≈ ΔErelax, could occur.19,20 Recently, it has been proposed that a modern solvation model, SMD, is sufficiently reliable that solution-phase Gibbs free energies can be obtained directly, using vibrational frequencies calculated in the solvent phase in conjunction with the partition functions derived for an ideal gas.20 It has also been proposed that the accurate calculations that are usually performed in the gas phase can give reliable results in combination with the SMD solvation model.21,22 In contrast, it is recommended that many other implicit solvation models only be applied in conjunction with the computational protocols used in their original parametrization.23 In the present work, we evaluate different approaches for conformational searching and the calculation of solution-phase free energies of flexible molecules. As Gibbs free energies of solvation are rarely measured in isolation, here we consider their use for the prediction of pKa values of mono- and polyprotic amines as a case study. pKa values are important in their own right for the insight they provide into the ability of species to engage in proton transfer reactions and their associated activity in different environments.24 pKa calculations for polyprotic amines are expected to be particularly challenging due to their multiple charges, conformational flexibility and the role of hydrogen bonding and other explicit solute−solvent interactions in determining their solvation energies. Herein we compare various protocols for conformational searching and solution-phase Gibbs free energy determination, and different approaches to the subsequent pKa prediction, with a view to making more general recommendations about the treatment of flexible and multiply charged molecules in solution.

Exploration of Conformational Space. All possible molecular conformations were constructed for rotations of 120° around each sp3−sp3 bond and 60° around each sp2−sp3 bond. 60° increments were also used for rotations around NH··· OH2 hydrogen bonds. Where the conformational space was very large (>81 structures), the energy-directed tree search algorithm (EDTS) was used to guide the search for the most stable conformation.26 Geometries and Frequencies. Molecular geometries were optimized and harmonic vibrational frequencies determined using density functional theory (see Supporting Information for full details). Where explicit solvation was included, a single water molecule was associated with each charged site. (I.e. no water molecules for neutral amines, one for monocations and two for dications.) No geometric constraints were applied, allowing bridged structures to form where appropriate. Geometries of the conformations identified as being the most stable for each species are reported in Table S10. Partition functions, and hence entropies and thermal corrections, were calculated using standard textbook formulas for the statistical thermodynamics of an ideal gas under the harmonic oscillator rigid rotor approximation at 298 K. Vibrational frequencies were scaled using recommended scaling factors for the calculation of zero point energies, thermal corrections and entropies.27 The value of using the quasi-harmonic oscillator approximation20 to minimize the effects of errors in very low frequency vibrations also was explored. Gibbs Free Energy Calculations. Four approaches to the calculation of solution-phase free energies were employed in this

■

METHODS Test Set. Accurate experimental base association constants are available for a wide range of amine-containing species.25 These data are an excellent resource for investigating the effects of molecular conformation choice on pKa predictions for the corresponding conjugate acids. A test set of 44 amine-containing species was constructed, comprising 20 primary amines, nine secondary amines, five tertiary amines and 10 species containing two amine groups. Species with long chains, branched chains and ring systems were included. See Chart 1 for structures and Table S1 for experimental pKas. 5218

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A work. These were based on two model chemistries: a comparatively inexpensive method (LL), based on M05-2X28/ cc-pVTZ29 energies, was used to predict approximate solutionphase Gibbs free energies (Gsoln * ) for initial comparison of conformer stabilities; while a higher-level (HL) approach, based on G3(MP2,CC,+),30 was used to calculate more accurate G*soln data for the prediction of pKas. In both cases, Gsoln * values could be determined using both entirely solution-phase calculations (LLsol and HL-sol), and via the use of a thermocycle (LL-TC and HL-TC). The SMD solvation model9 was used for all solutionphase calculations. Full details of the model chemistries are provided in the Supporting Information. Gsoln * values calculated via the LL-sol, HL-sol, and HL-TC approaches, along with their components, are reported in Tables S7, S8, and S9, respectively. pKa Calculation Schemes. Much study has gone into identifying reaction schemes which lead to accurate pKa predictions; see, for example, refs 15 and 31−37. Among the many schemes reported, three broad classes can be identified: direct methods, methods involving the donation of a proton to a reference acid, and methods in which charged sites are clustered with one or more explicit solvent molecules. It is not the goal of this work to reassess the merits of these various approaches, but rather to highlight additional concerns associated with the use of these schemes to investigate pKas of large, flexible, and multiply ionizable species. As such, we have selected four representative schemes: one example of each of the three classes, plus a combination method involving the use of both a reference species and explicit solvation of both the reference acid and species of interest. This final approach has previously been shown to lead to superior accuracy in pKa calculations.31 We label these the direct, proton exchange (PEX), implicit−explicit (IE) and IE+PEX methods. Full details how to determine pKa via each of these approaches are presented in the Supporting Information. pKas calculated using the LL-sol, HL-sol, and HL-TC approaches with each of the four pKa calculation schemes are reported in Table S5. Deviations of these values from experiment can be found in Table S6.

Figure 1. Relative stabilities of selected low lying conformations of 1,3propanediamine determined using the Egas, Esoln, Ggas, and Gsoln metrics and the LL model chemistry (here the LL-sol approach is used for Gsoln). Structures for four important conformations are shown, with the color of the box corresponding to the trace in the plot.

small, it would still contribute an error of 0.3 pKa units to the final pKa prediction. Clearly, gas-phase metrics are not appropriate for conformational searching, and comparison of solution-phase Gibbs free energies for low lying conformations is highly recommended. This raises the question of how Gsoln should be calculated. Until recently, thermocycle based methods were accepted to be the most reliable approach, or at least the approach most consistent with the original parametrization of the contributing calculations.19 As described in the Supporting Information, the gasphase free energy of each conformer is corrected by the free energy released via solvation (ΔGsolv); the free energy gained due to structural relaxation in the solvent can also be (approximately) accounted for via the corresponding energy change (ΔErelax). In Figure 2, the contributions of ΔGsolv and ΔErelax to Gsolv(LL-TC) are shown for the conformations of the 1,3-propanediamine dication relative to their values for the fully extended species (red). Absolute values for these corrections can be found in

■

RESULTS AND DISCUSSION Choosing Conformations. The usefulness of the gas-phase energy, Egas, solution-phase energy, Esoln, and gas-phase Gibbs free energy, Ggas, as metrics for identifying the most stable conformation of a molecule in solution (lowest Gsoln) was assessed using 1,3-propanediamine as a test case. Geometries for all possible molecular conformations were optimized with M052X/6-31G(d). Electronic energies and free energies were calculated using the LL model chemistry. The relative stabilities of selected conformations are shown in Figure 1; corresponding values are reported in Table S2. For each metric, the (Gibbs free) energy of each conformation is shown relative to the conformation which is predicted to be most stable. In all cases, at least the three most stable conformations are shown. As is clear from the figure, conformations which are predicted to be stable in the gas phase are high in (Gibbs free) energy in solution. While the conformation shown in dark purple is the most stable according to the Egas and Ggas metrics, it is actually one of the highest free energy conformations in the Gsoln manifold (+10.8 kJ mol−1). Although this conformation has an intramolecular hydrogen bond, the conformation shown in light purple, which has a similar profile, does not. Likewise, the most stable conformation according to the Esoln metric (green) lies 1.7 kJ mol −1 above the lowest G soln conformation. Although this free energy difference is relatively

Figure 2. Relative contributions of the components of Gsoln,TC to the final value for different conformations of the 1,3-propanediamine dication. Values are reported relative to those of the red conformation (the conformation with the smallest magnitude corrections) in all cases; corresponding absolute values can be found in Table S3. Gsoln,soln values for the four conformations are provided for comparison. Solution-phase structures of each conformation are shown in CPK coloring; these are superimposed on the corresponding gas-phase structures, shown in color. (Note: the purple solution-phase structure relaxes to the blue structure in the gas phase.) All calculations performed using the LL model chemistry. 5219

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A

Given these concerns in relation to thermocycle-based methods and the claimed reliability of modern solvation models, the possibility of calculating Gsoln using an entirely solution-phase methodology is appealing. Unfortunately, this approach also has its weaknesses. This is highlighted by the cation of 1,3propanediamine, for which both Gsoln(LL-sol) and Gsoln(HLsol) predict a conformation with a tight 6-membered ring, formed by an intramolecular hydrogen bond, to be the most stable (by 4.0 and 2.5 kJ mol−1, respectively; see Figure S2). However, when this conformation is used to predict the pKas for the first and second protonation of 1,3-propanediamine, the first is overpredicted by 0.9 pKa units and the second underpredicted by 0.7 units. When pKas are calculated using the most stable 1,3propanediamine conformer that has no intramolecular hydrogen bond, both results move closer to the experimental values (now +0.5 and −0.2 pKa units, respectively). This trend is observed for all diamine species where an intramolecular hydrogen bond is predicted, indicating that the solution-phase methodology is significantly overstabilizing these conformations with respect to their non-hydrogen-bonded counterparts. In general, we find that the Gsoln,soln approach provides the most reliable metric for identifying appropriate conformations for pKa calculations. Care must be taken, however, when there is the potential for formation of tight intramolecularly hydrogenbonded loops. In these cases, comparison of multiple approaches, and ideally benchmarking using related systems, is required. We have also assessed whether using a Boltzmann average of the free energies of all conformations improves pKa predictions; this was found to give no improvement to the results. Throughout the rest of this work, low lying conformers of each species were initially identified based on their Gsoln values calculated according to the LL-sol methodology. The HL-sol approach was then applied to each of these structures to identify the true lowest Gsoln conformation. Where the conformation involved a tight intramolecularly hydrogen-bonded loop, the lowest Gsoln(HL-sol) conformation with no intramolecular hydrogen bonding was used. Dealing with Frequencies. In order to calculate Gibbs free energies, the vibrational frequencies of the molecule are required. While the zero-point energy of a molecule is most sensitive to its highest frequency vibrations, the thermal corrections to the enthalpy and the vibrational entropy of the molecule depend most heavily on the lowest energy modes. Unfortunately, these low energy modes become more prevalent as the molecular size increases. Additionally, they are notoriously difficult to calculate accurately, and even small errors in very low energy modes can lead to large inaccuracies in free energies and hence pKas. In order to reduce the impact of these errors, the quasi-harmonic oscillator (QHO) approach was proposed.20 This assumes that low energy modes should be relatively unaffected by protonation, and should therefore not contribute to the free energy of reaction. As a result, low energy frequencies can be systematically scaled up a consistent value; in the original formulation, all frequencies below 100 cm−1 were scaled to 100 cm−1. For our test set, we used a very rigorous procedure to ensure maximal accuracy in our prediction of low vibrational frequencies. This allowed the impact of the QHO approach to be assessed and the most reliable cutoff value in various situations to be identified. When using the direct and PEX pKa calculation schemes (no explicit H2O present), application of the QHO approach with varying cutoff values did not generally have a significant effect on the pKa values predicted for our test set (see Figure S3). This suggests that any cutoff value up to 100 cm−1 is appropriate.

Table S3. Clearly, the relative contributions of the both the solvation and relaxations terms vary greatly between conformations, each by up to 50 kJ mol−1 in this case. In fact, the differences in the correction terms between conformations swamp the variations in the free energies themselves. This raises concerns that small percentage errors in either of these terms could change the prediction of the most stable conformer. In particular, the approximation that ΔGrelax ≈ ΔErelax is called into question. As can be seen in Figure 2, both the blue and green conformations show appreciable differences between their gasand solution-phase geometries. This is due to the two charged groups being partially shielded from each other by the dielectric continuum in the solution-phase calculations, but experiencing increased electrostatic repulsion in the gas phase. Associated with these geometric changes will be changes in the vibrational frequencies, and hence the vibrational contributions to the free energies. Unfortunately, it is not possible to estimate the magnitude or even the net sign of these effects. However, chemical intuition would suggest that it would be surprising if the green and blue conformations, with two NH3+ groups separated by 4.1 and 4.4 Å, respectively, were truly more stable in aqueous solution than the red conformation, for which the two charged groups are at maximal separation (4.9 Å). Additional complications arise for thermocycle based methods because not all conformations that are stable in solution are also stable in the gas phase. For example, the purple structure in Figure 2 is identified as a distinct conformation in aqueous solution, however upon optimization in the gas phase the electrostatic repulsion between the two NH3+ groups forces a conformational change to the blue conformer. Hence, the gasand solution-phase calculations of the thermocycle actually involve different conformations. While this conformer is not predicted to have the lowest Gsoln for this system, it is possible that in other cases structures undergoing analogous conformational changes may be predicted to be the most stable. Large structural changes between phases are also seen when the system involves weak interactions, e.g. in transition states,18 or where hydrogen bonded interactions with explicit solvent molecules are present. For protonated amine−water complexes, the plane of the H2O molecule sits at ∼155° to the hydrogen bond in the gas phase, while in solution this angle is closer to 115° (see Figure S1). pKa predictions provide a useful tool to assess whether thermocycle-based methods can reliably determine Gsoln when there are large structural differences between the gas- and solution-phase geometries. As can be seen in Table S4 (and also Figures S4c and S5), even when using the HL model chemistry, thermocycle methods perform very poorly for pKa predictions where the large structural changes associated with RNH+···H2O complexes are involved. In particular, the IE HL-TC approach gives by far the worst results of all approaches described in this work, with a mean absolute deviation from experiment of 0.9 pKa units and a maximum error of 4.3. Hence, thermocycle methods can be unreliable when the system undergoes significant structural or conformational changes between the gas and solution phases. This is a particular concern when performing conformational searches, for which accurate relative free energies of all conformations are needed. The computational cost of calculating Gsoln,TC (involving two geometry optimizations plus the determination of ΔGsolv) also makes the use of thermocyclebased approaches unattractive for systematic conformational searches. 5220

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A

Figure 3. Effect of varying the QHO cutoff value on pKa predictions for the protonation of primary amine sites of species in the amine test set. Deviations from experiment of pKas calculated using the IE reaction scheme with the HL-sol approach are shown in pKa units. Red gridlines indicating deviations of ±1 unit are included to guide the eye. Mean signed errors (MSE), mean unsigned errors (MUE), standard deviations (SD) and maximum errors (MaxE) in pKa units are shown in each case. While varying the QHO cutoff values between 20 and 50 cm−1 (magenta, green, yellow and blue) does not significantly affect either the pKa values predicted or their overall accuracy, a cutoff of 100 cm−1 (red) generally produces an appreciable underestimate of the experimental pKa. This is particularly observed in the doubling of the MUE.

used in this work, ammonia (P1) or methylamine (P2) would be common choices. As shown in Figure 4, however, using these species as the references for predicting the pKas of large flexible amines can lead to significant errors. Over the entire test set, using the direct method (green) gives a mean deviation from experimental of −0.13 pKa units. Using the proton exchange method with methylamine as the reference (PEX(P2), yellow) produces a vertical shift in the pKa predictions; this actually increases the mean deviation from experiment to +0.38. If ammonia was used as the reference, the vertical shift would be even more dramatic, and the agreement with experiment even poorer. In contrast, using a medium sized amine, such as pentylamine (P9), as the reference (PEX(P9), blue) reduces the mean deviation from experiment to +0.08 units. Hence, using a reference acid of similar size and flexibility to the acid of interest gives better pKa predictions than simply using the smallest related species. (It is important to note that PEX schemes yield relatively small improvements over the direct method for this test set simply because the direct method itself gives a relatively good description of amines. For other classes of molecules, the improvements can be much more significant. For instance, for anilines the direct method generally predicts pKas that are ∼2 units too low.) Further examination of Figure 4 reveals that many of the systems that still show errors of more than one pKa unit involve the protonation of secondary or tertiary amine sites. Refinements could therefore be made by using different reference species for primary, secondary and tertiary amines. If small amines of the appropriate class are used, e.g. methylamine (P2) for primary amines, dimethylamine (S1) for secondary amines and trimethylamine (T1) for tertiary species (PEX(P2/S1/T1), purple), the standard deviation from experiment across the test set falls to 0.41 pKa units, in comparison with 0.54 for the direct method, 0.63 for PEX(P2) and 0.54 for PEX(P9). Choosing a larger, more

When explicit water molecules are included in the calculation (IE and IE+PEX schemes), however, the cutoff value becomes much more important. pKa predictions calculated using the IE method and various QHO cutoffs are shown in Figure 3. For some systems, the pKa is seen to vary by nearly two pKa units solely due to the cutoff value used (e.g., P15, octylamine). This is an excellent example of why a good choice of the QHO cutoff is required. The appreciable shift between cutoffs of 0 and 25 cm−1 indicates that at least one low energy mode has significantly different values in the neutral and protonated species. The use of a 100 cm−1 QHO cutoff, however, leads to a systematic underestimation of the pKa for almost all species. This is because the NH···OH2 bending modes, which are present in the hydrated cation but not in the neutral species, are being neglected due to the scaling. Comparing the average statistics indicates that cutoffs of 25, 30, and 40 cm−1 all give comparable accuracy for pKa prediction; we have thus chosen a value of 30 cm−1 for species both with and without explicit solvation throughout the rest of this work. In smaller scale investigations, we recommend examining the low energy modes individually and only scaling those which should not be affected by the reaction in question. Proton Exchange Schemes: Choosing a Good Reference Acid. Proton exchange schemes attempt to minimize systematic errors in pKa predictions. In essence, they involve assessing the error in a direct method pKa calculation for an acid with an accurately known pKa, then using this to correct the pKa prediction for the acid of interest. It is well-known that the sign and magnitude of the correction depends strongly on the class of acid being considered, and that it is necessary to choose a reference acid that is “structurally similar” to the acid being studied.31 Often the smallest member of a particular class of molecules is chosen to be the reference;31,32 for instance, for the amine test set 5221

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A

Figure 4. Effect of varying the reference species on pKa predictions for the amine test set. Deviations from experiment of pKas calculated using the PEX reaction scheme with the HL-TC approach are shown in pKa units. Red gridlines indicating deviations of ±1 unit are included to guide the eye. Mean signed errors (MSE), mean unsigned errors (MUE), standard deviations (SD), and maximum errors (MaxE) in pKa units are shown in each case.

flexible species of the relevant class for the reference acid gives further improvement. If pentylamine (P9) is used for primary amines, diethylamine (S2) for secondary amines and diethylmethylamine (T3) for tertiary amines (PEX(P9/S2/T3), red), the standard deviation falls to 0.38. Significant improvement is also observed in the mean absolute deviation (MUE): 0.38 for PEX(P9/S2/T3) vs 0.59 for the direct method, 0.81 for PEX(P2), 0.56 for PEX(P9), and 0.52 for PEX(P2/S1/T1). With the PEX(P9/S2/T3) approach, only two examples are seen where the absolute error is significantly above 1 pKa unit: ammonia (P1), with Td symmetry and no methyl group; and pyrrolidine (S5), with a strained five-membered ring. In both cases, it could be argued that the species is sufficiently different from the rest of the class that the chosen reference may not be appropriate. Importantly, the corrections obtained for primary, secondary and tertiary monoamines apply equally well to the corresponding amine sites in polyamines. In summary, when choosing a reference species for a complex acid it is important to consider not only the chemistry of the species, but to also seek, where possible, to match the size, flexibility and structural strain of the molecule. Unless otherwise noted, where the PEX and IE-PEX schemes are used throughout this work, the combination of P9, S2, and T3 reference species is used. Comparison of pKa Calculation Approaches. A comparison of the performance of the various pKa calculation approaches considered in this work for the amine test set is provided in Figures S4 and S5. Overall statistics can be found in Table S4. The best performing approaches for various scenarios are shown in Figure 5. Of particular interest with regard to the calculation of pKas for large, flexible molecules is the performance of the less computationally expensive LL model chemistry. For most species, the LL-sol approach is found to give similar pKa predictions to the corresponding HL-sol methodology. Inclusion of an explicit water molecule (the IE reaction scheme) further

improves the results of the LL-sol approach, such that it actually out-performs HL-sol for tertiary amine species. In fact, of all the approaches that can be used when no reference data are available (direct and IE schemes), the IE LL-sol methodology gives the best overall performance for the monoamines in this test set (mean error −0.03 pKa units; standard deviation 0.21; mean absolute deviation 0.39; max error 1.6). It is important to bear in mind, however, that this good agreement may be due to the primary monoamines being a particularly easy case for pKa calculation, and such good accuracy may not be obtained for other classes of molecules. Indeed, the LL-sol methodology did struggle to reliably describe diamines, the specific interest of this paper. In particular, poor pKa predictions (error >1 pKa unit) are usually obtained for the second protonation step, most significantly when the two protonation sites are in close proximity (D1, D2, and D6−D10). Indeed, the HL-sol approach also finds the second protonation of diamines challenging unless an explicit water molecule is included in the calculation. On the basis of the results from this test set, the most reliable non-PEX approach for species with two ionizable groups in close proximity is the direct HL-TC methodology. The inclusion of an explicit water molecules in pKa calculations (IE or IE+PEX schemes) gave variable responses. For the LL-sol model chemistry, errors in the pKa predictions were reduced for monoamines, but not for diamines. In contrast, when the IE HLsol approach is used, improvements are seen in the second protonation of diamines (due to a better description of the singly charged species leading to reduced negative deviations from experiment), however the agreement with experimental pKas becomes worse for the first protonation step, and also for triamine species. It is only when a proton exchange scheme is also applied (IE+PEX HL-sol) that good overall agreement with experiment is obtained. (Note: It is important to remember that HL-TC methodology cannot be applied when explicit water molecules are present, as differences in gas- and solution-phase 5222

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A

Figure 5. Deviations from experimental pKas (in pKa units) for the best performing approaches under different circumstances. Blue: IE LL-sol approach for pKas of monoamines where no reference data are available. Green: Direct HL-TC approach for pKas of diamines where no reference data are available. Yellow: PEX HL-TC approach for all species where inclusion of explicit water molecules is not necessary/desired. Red: IE+PEX HL-sol approach for all species where inclusion of explicit water molecules is required. Red gridlines indicating deviations of ±1 unit are included to guide the eye. Mean signed errors (MSE), mean unsigned errors (MUE), standard deviations (SD), and maximum errors (MaxE) in pKa units are shown in each case.

• It is important that the conformation with the lowest solution-phase free energy (Gsoln) be identified. • Thermocycle-based approaches to calculating Gsoln are not recommended for this purpose. • Care must be taken for systems where the conformation predicted to be most stable shows strong intramolecular interactions. Concerning frequencies: • It is important to take steps to avoid large errors due to inaccuracies in the prediction of low energy vibrational modes. • The quasi-harmonic oscillator approximation is a reliable approach for dealing with this issue. • Cutoff values between 20 and 50 cm−1 lead to good pKa predictions for this test set. • Where possible, modes should be examined individually, and only those not expected to contribute to the reaction should be scaled. Concerning reference acids for proton exchange (PEX) schemes: • Where possible, the chemical environment surounding the ionizable site should be considered; i.e., use different reference species for primary, secondary, tertiary and strained amine sites. • The size of the reference acid is also relevant. The smallest species of a particular class can often carry different errors to larger molecules in the same series. Concerning model chemistries and reaction schemes: • Approaches where only density functional theory is employed (e.g., LL-sol) can give good pKas, so long as the chemistry around the ionization site is not too complex. When explicit water molecules were included,

geometries (see Figure S2 and the Discussion) can lead to large errors in the pKa predictions.) As expected, the results from this test set show that, where possible, the use of a PEX scheme is to be preferred. The two methodologies which provided the best agreement with experiment across the entire amine test set were PEX HL-TC and IE+PEX HL-sol. (PEX HL-TC: mean error −0.10 pKa units; standard deviation 0.38; mean absolute deviation 0.21; max error 1.6 for pyrrolidine, as described earlier. IE+PEX HL-sol: mean error −0.11 pKa units; standard deviation 0.36; mean absolute deviation 0.21; max error 1.2, again for pyrrolidine.) On the basis of the results from this test set, for systems where the addition of explicit water molecules is not feasible or desired, the PEX HL-TC approach is recommended. Where explicit solvent is required to model the chemistry of the system, however, use of the IE+PEX HL-sol methodology is advised. For systems with only one functional group, the IE LL-sol approach was found to generally give acceptable results without the use of a reference species. For cases with multiple ionizable sites in close proximity, the best non-PEX methodology tested here was the direct HL-TC approach. It is important to recognize that the comparative performance of different methodologies found for this test set may not translate to other classes of molecules.

■

SUMMARY Extensive work by many different research groups over the past 15 years has provided a good understanding of the factors which must be addressed when calculating pKas of small, rigid molecules. Here we extend this work to identify additional concerns which must be considered for large, flexible systems, including those containing multiple ionizable groups. We find that Concerning conformations: 5223

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A Notes

this approach performed remarkably well for monoamines, giving the best non-PEX pKa predictions for these species. • Thermocycle-based approaches give reliable pKas for calculations involving bare ions (particularly PEX HLTC), however they should not be applied when explicit solvent molecules are present. • Performing high-level calculations in the solution phase (HL-sol) gives comparable results to HL-TC when a PEX reaction scheme is used, however this approach has the advantage of also being able to be applied in conjunction with explicit solvation. • Overall, including explicit solvation via the IE+PEX HLsol methodology gave the most reliable pKa predictions for the molecules in this test set. Concerning explicit and implicit solvation: • Most species in the test set were well described without introducing explicit water molecules; it was only in particularly challenging cases that the added cost of including explicit solvation yielded improved results, and then only for IE+PEX HL-sol. • On the basis of our results, it appears that explicit solvation should not be included when a thermocycle method is employed. • Note: Combining high-level or large basis set calculations with implicit solvation is only recommended for the SMD solvation model. The authors of other solvation models, such as PCM and COSMO-RS, have previously recommended that only the computational protocol used in the parametrization of the model be employed for best results.23

■

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Generous allocations of supercomputing time on the National Facility of the Australian National Computational Infrastructure are gratefully acknowledged.

■

(1) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999−3094. (2) Miertuš, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilizaion of ab initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117−129. (3) Foresman, J. B.; Keith, T. A.; Wiberg, K. B.; Snoonian, J.; Frisch, M. J. Solvent effects. 5. Influence of cavity shape, truncation of electrostatics, and electron correlation on ab initio reaction field calculations. J. Phys. Chem. 1996, 100, 16098−16104. (4) Pascual-ahuir, J. L.; Silla, E.; Tuñon, I. GEPOL: An improved description of molecular surfaces. III. A new algorithm for the computation of a solvent-excluding surface. J. Comput. Chem. 1994, 15, 1127−1138. (5) Miertus̃, S.; Tomasi, J. Approximate evaluations of the electrostatic free energy and internal energy changes in solution processes. Chem. Phys. 1982, 65, 239−245. (6) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model. J. Comput. Chem. 2003, 24, 669−681. (7) Klamt, A.; Schuurmann, 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. (8) Klamt, A. Conductor-Like Screening Model for Real Solvents - a new approach to the quantitative calculation of solvation phenomena. J. Phys. Chem. 1995, 99, 2224−2235. (9) 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. (10) Soteras, I.; Curutchet, C.; Bidon-Chanal, A.; Orozco, M.; Luque, F. J. Extension of the MST model to the IEF formalism: HF and B3LYP parametrizations. J. Mol. Struct.: THEOCHEM 2005, 727, 29−40. (11) Mobley, D. L.; Wymer, K. L.; Lim, N. M.; Guthrie, J. P. Blind prediction of solvation free energies from the SAMPL4 challenge. J. Comput.-Aided Mol. Des. 2014, 28, 135−150. (12) Gupta, M.; da Silva, E. F.; Svendsen, H. F. Modeling temperature dependency of amine basicity using PCM and SM8T implicit solvation models. J. Phys. Chem. B 2012, 116, 1865−1875. (13) Sutton, C. C. R.; Franks, G. V.; da Silva, G. First principles pKa calculations on carboxylic acids using the SMD solvation model: effect of thermodynamic cycle, model chemistry, and explicit solvent molecules. J. Phys. Chem. B 2012, 116, 11999−12006. (14) Chen, Y.-L.; Barlow, D. J.; Kong, X.-L.; Ma, Y.-M.; Hider, R. C. Prediction of 3-hydroxypyridin-4-one (HPO) hydroxyl pKa values. Dalton Trans. 2012, 41, 6549−6557. (15) Sastre, S.; Casasnovas, R.; Munoz, F.; Frau, J. Isodesmic reaction for accurate theoretical pKa calculations of amino acids and peptides. Phys. Chem. Chem. Phys. 2016, 18, 11202−11212. (16) Zhang, S.; Baker, J.; Pulay, P. A reliable and efficient first principles-based method for predicting pKa values. 2. Organic acids. J. Phys. Chem. A 2010, 114, 432−442. (17) Zhang, S. A reliable and efficient first principles-based method for predicting pKa values. III. Adding explicit water molecules: Can the theoretical slope be reproduced and pKa values predicted more accurately? J. Comput. Chem. 2012, 33, 517−526. (18) Lin, C. Y.; Izgorodina, E. I.; Coote, M. L. First principles prediction of the propagation rate coefficients of acrylic and vinyl esters: are we there yet? Macromolecules 2010, 43, 553−560.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.7b04133. Full details of the methodologies used in this work, molecular structures that highlight differences in geometry predictions obtained from different computational approaches, molecular geometries, energies, and free energies (including components thereof) calculated via the different approaches described in this work, and pKas of the molecules in the test set derived from experiments and calculations, including deviations from experiment and statistics assessing the performance of the different computational approaches (PDF).

■

REFERENCES

AUTHOR INFORMATION

Corresponding Authors

*(N.L.H.) E-mail: [email protected]. Telephone: (+61) 2 6125 4039. *(M.L.C.) E-mail: [email protected]. Telephone: (+61) 2 6125 3771. ORCID

Naomi L. Haworth: 0000-0002-3299-3137 Michelle L. Coote: 0000-0003-0828-7053 Funding

This work was funded via the Australian Research Council Centre of Excellence for Electromaterials Science (CE140100012). 5224

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225

Article

The Journal of Physical Chemistry A (19) Ho, J.; Klamt, A.; Coote, M. L. Comment on the correct use of continuum solvent models. J. Phys. Chem. A 2010, 114, 13442−13444. (20) 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. (21) Casasnovas, R.; Fernández, D.; Ortega-Castro, J.; Frau, J.; Donoso, J.; Muñoz, F. Avoiding gas-phase calculations in theoretical pKa predictions. Theor. Chem. Acc. 2011, 130, 1−13. (22) Ho, J. Are thermodynamic cycles necessary for continuum solvent calculation of pKas and reduction potentials? Phys. Chem. Chem. Phys. 2015, 17, 2859−2868. (23) Klamt, A.; Mennucci, B.; Tomasi, J.; Barone, V.; Curutchet, C.; Orozco, M.; Luque, F. J. On the performance of continuum solvation methods. A comment on “Universal approaches to solvation modeling. Acc. Chem. Res. 2009, 42, 489−492. (24) Ho, J.; Coote, M. L. First-principles prediction of acidities in the gas and solution phase. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 649−660. (25) Haynes, W. M.; Lide, D. R.; Bruno, T. J. CRC Handbook of Chemistry and Physics: a Ready-reference Book of Chemical and Physical Data; CRC Press: Boca Raton, FL, 2015. (26) Izgorodina, E. I.; Yeh Lin, C.; Coote, M. L. Energy-directed tree search: an efficient systematic algorithm for finding the lowest energy conformation of molecules. Phys. Chem. Chem. Phys. 2007, 9, 2507− 2516. (27) Merrick, J. P.; Moran, D.; Radom, L. An evaluation of harmonic vibrational frequency scale factors. J. Phys. Chem. A 2007, 111, 11683− 11700. (28) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of density functionals by combining the method of constraint satisfaction with parametrization for thermochemistry, thermochemical kinetics, and noncovalent interactions. J. Chem. Theory Comput. 2006, 2, 364−382. (29) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (30) Henry, D. J.; Sullivan, M. B.; Radom, L. G3-RAD and G3X-RAD: Modified Gaussian-3 (G3) and Gaussian-3X (G3X) procedures for radical thermochemistry. J. Chem. Phys. 2003, 118, 4849−4860. (31) Ho, J.; Coote, M. A universal approach for continuum solvent pKa calculations: are we there yet? Theor. Chem. Acc. 2010, 125, 3−21. (32) Rebollar-Zepeda, A. M.; Galano, A. First principles calculations of pKa values of amines in aqueous solution: Application to neurotransmitters. Int. J. Quantum Chem. 2012, 112, 3449−3460. (33) Thapa, B.; Schlegel, H. B. Density functional theory calculation of pKas of thiols in aqueous solution using explicit water molecules and the polarizable continuum model. J. Phys. Chem. A 2016, 120, 5726−5735. (34) Eckert, F.; Diedenhofen, M.; Klamt, A. Towards a first principles prediction of pKa: COSMO-RS and the cluster-continuum approach. Mol. Phys. 2010, 108, 229−241. (35) Pliego, J. R.; Riveros, J. M. Theoretical Calculation of pKa using the cluster−continuum model. J. Phys. Chem. A 2002, 106, 7434−7439. (36) 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. (37) 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.

5225

DOI: 10.1021/acs.jpca.7b04133 J. Phys. Chem. A 2017, 121, 5217−5225