Dissociation constants of weak acids from ab initio molecular dynamics using metadynamics: influence of the inductive effect and hydrogen bonding on pKa values.

Subscriber access provided by UNIV OF TOLEDO

Article

Dissociation Constants of Weak Acids from Ab Initio Molecular Dynamics Using Metadynamics: the Influence of the Inductive Effect and Hydrogen Bonding on pKa Values Sukumaran Vasudevan, and Anil Kumar Tummanapelli J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp5088898 • Publication Date (Web): 06 Nov 2014 Downloaded from http://pubs.acs.org on November 9, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Dissociation Constants of Weak Acids from Ab Initio Molecular Dynamics Using Metadynamics: the Influence of the Inductive Effect and Hydrogen Bonding on pKa Values

Anil Kumar Tummanapelli and Sukumaran Vasudevan* Department of Inorganic and Physical Chemistry Indian Institute of Science, Bangalore 560012, INDIA

* Author to whom correspondence may be addressed. E-mail: [email protected]. Tel: +91-80-2293-2661. Fax: +91-80-2360-1552/0683;

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 25

Abstract The theoretical estimation of the dissociation constant, or pKa, of weak acids continues to be a challenging field. Here we show that ab initio Car-Parrinello molecular dynamics simulations in conjunctionwith metadynamics calculations of the free energy profile of the dissociation reaction provide reasonable estimates of the pKa value. Water molecules, sufficient to complete the three hydration shells surrounding the acid molecule, were included explicitly in the computation procedure. The free-energy profiles exhibit two distinct minima corresponding to the dissociated and neutral states of the acid and the difference in their values provides the estimate for pKa. We show for a series of organic acids that CPMD simulations in conjunction with metadynamics can provide reasonable estimates of pKa values. The acids investigated were aliphatic carboxylic acids, chlorine substituted carboxylic acids, cis- and transbutenedioic, and the isomers of hydroxy-benzoic acid.

These systems were chosen to highlight that the procedure could

correctly account for the influence of the inductive effect as well as hydrogen bonding on pKa values of weak organic acids. In both situations the CPMD-metadynamics procedure faithfully reproduces the experimentally observed trend and the magnitudes of the pKa values.

Keywords: Weak acid dissociation, Car-Parrinello Molecular Dynamics, Metadynamics, Free Energy calculations.


2

Page 3 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION The pKa value, or acid dissociation constant, is a fundamental property of weak acids that holds the key to understanding their reactivity. Its experimental determination, at least for monoprotic acids or in polyprotic acids where successive pKa values are well separated, is easily realized by a relatively straight-forward application of the Henderson-Hasselbalch equation, and this exercise

often forms a part of most undergraduate laboratory curricula. The theoretical or

computational estimation of pKa values, even for the simplest acids, is

however, a non-trivial

exercise.1-3 The need for a theoretical estimate of pKa may not be apparent for simple monoprotic acids but in situations where there are multiple ionizable groups or in the case of large systems, like proteins, where ionizable residues are buried deep within, the ability to predict pKa values is clearly desirable. Here we show that ab initio Car-Parrinello Molecular Dynamics (CPMD)4 in conjunction with metadynamics computation5-8 of the free-energy landscape provides a simple method to estimate pKa values of weak organic acids with reasonable accuracy. We demonstrate, using text book examples, that the procedure correctly predicts the effect of the inductive effect as well as hydrogen bonding, two well known factors that affect acid dissociation, and show for a series of substituted and un-substituted aliphatic and aromatic carboxylic acids that the method provides estimates of the pKa values in reasonable agreement with experimental data. Current approaches for theoretically estimating pKa values usually involve constructing the appropriate thermo-chemical acid dissociation cycle followed by evaluation of the free-energy change at each step of the cycle.1 The calculations are demanding, for an error of 1 kcal/mole in the free energy value results in a error of 0.74 pKa units.9 While changes in free-energy for the steps in the gas-phase can be evaluated with sufficient accuracy the same is not true when solvation effects have to be considered.10 By far the most challenging aspect in the theoretical


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 25

prediction of pKa values is to account for the solvent in a physically meaningful way. This is both conceptually and computationally demanding. The majority of these calculations have been carried out by treating the solute at the quantum mechanical (QM) level and

modeling the

solvent as a dielectric continuum.11 The drawback of these procedures, and a major source of error, is that short-range intermolecular interactions with solvent molecules, such as hydrogen bonds and ion-dipole interactions, are not explicitly accounted for although they are particularly important as the calculations involve ionic species. A number of different strategies have been developed to improve the accuracy of the calculations by explicit inclusion of solvent water molecules in the acid dissociation process.12-14 For strongly coupled solute-solvent systems explicit inclusion of water molecules in the quantum mechanical part of the calculation can greatly improves the accuracy in computed pKa values, especially for certain continuum models, such as the SM6 solvation model.3 A number of hybrid approaches have been developed to capture the first and second solvation directed effects in the context of continuum models. These have been summarized in recent reviews.3,15-18 An alternate approach to the continuum models is to model the solvent molecules using molecular mechanics (MM) using an appropriate force fields while describing the bond-breaking by QM procedures.19-21 Instead of constructing thermo-chemical cycles, free energies are obtained by constrained QM-MM molecular dynamics in conjunction with umbrella sampling.22 These methods are able to predict pKa values of small molecules, e.g., HF, HCOOH, CH3COOH and H2CO3 in reasonable agreement with experiment.22 It has been shown that boundary issues inherent to the QM-MM method can be avoided using reactive-MD that are appropriately parameterized. Multiscale reactive–MD simulations have been used to predict absolute pKa values of aspartic and glutamic amino acids.23


4

Page 5 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A key issue in any dissociation event is how the solvating water molecules arrange themselves both spatially and dynamically around the ionizable and ionized species. This, in principle, is best described by ab initio or quantum molecular dynamics molecules in the ensemble are treated explicitly.8

wherein all solvent

In these simulations the electronic structure

calculations are usually based on density functional methods and although these are more restricted as compared to the methods used in the hybrid continuum models have the advantage that solvation effects are handled with greater rigor. These methods are computationally more demanding but with increased availability of computer resources this is no longer a major constraint and consequently a number of investigations of acid dissociation using ab initio molecular dynamics including Car-Parrinello molecular dynamic have been reported.24-27 The systems investigated were HCOOH, CH3COOH and H2CO3.24-27 These simulations do not directly provide pKa values but procedures have been developed to obtain the acid dissociations constants. Thermodynamic free energies associated with the dissociation reaction have been calculated using the distributions of the vertical energy gap for the insertion or deletion of a proton.28 Using this approach, the pKa of various compounds, such as HCl, HCOOH, H2S and ortho and para salicylic acids have been calculated.28 The dissociation of histidine in aqueous media

has

been explored using constrained ab initio Car-Parrinello molecular dynamics

(CPMD).27 Free energies were obtained by integrating the potentials of mean force along a chosen coordinate for the deprotonation of histidine. Although the chosen reaction coordinate limited their analysis to stopping just beyond the transition state, the correct relative pKa value was obtained by subtracting the equivalent profile for a reference reaction, the auto-dissociation of water.


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 25

One of the factors that has inhibited the wide-spread use of ab initio MD methods to study the dissociation reaction, apart from demands on computational resources, is that dissociation of weak acids are rare events that require extremely long simulation times before one is observed, which further adds to computational costs. The metadynamics formalism provides a solution to this conundrum by preventing the system from revisiting regions of configuration space where it has been in the past.5 The formalism allows the system to escape the free-energy minima by biasing the dynamics with a history dependent potential (or force) that acts on select degrees of freedom, referred to as collective variables. The bias potentials, modeled by repulsive inverted Gaussians that are dropped during propagation, drive the system out of any free energy minima and allow it to explore the configurational space by a relatively quick and efficient sampling. CPMD in combination with metadynamics has been successfully used to sample the configurational space of acetic acid in bulk water.24 The dissociation reaction occurred along a well-characterized pathway and free energy minima corresponding to the associated and dissociated states of acetic acid were observed along with a shallow minimum in the free energy surface corresponding to a contact ion pair that is separated from the protonated state by a small energy barrier. Here we address the question as to whether it is possible to correctly estimate pKa values from the difference in the values of the free energy associated with the minima corresponding to the neutral and dissociated states of the acid in the metadynamics computed free-energy profile. Considering the nature of approximations both in the CPMD simulations and the metadynamics sampling procedure this may appear rather far-fetched. But, as we show here, since pKa values are obtained as the difference in free energy values, error magnitudes are reduced by cancellation. The focus of the present exercise was to investigate whether the CPMD-metadynamics formalism can not only provide reasonable estimates of the


6

Page 7 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dissociation constants but also capture the influence of the inductive effect and hydrogen bonding on the pKa values. We show here that the method predicts the pKa value of acetic acid in good agreement with experiment and then go on to show for a series of organic acids that ab initio CPMD simulations in conjunction with metadynamics calculations of the free energy profile of the dissociation reaction can provide reasonable estimates of the pKa value. The systems investigated were chosen to highlight whether the procedure could correctly account for the influence of the inductive effect as well as hydrogen bonding on pKa values of weak organic acids. SIMULATION METHODOLOGY The CPMD v3.11.2 program29 was used to perform the ab initio molecular dynamics simulation. Simulations were performed on an IBM Blue Gene cluster using 256 (128 × 2) processors. In the simulations, the organic acid were solvated by water molecules in a cubic box with periodic boundary conditions. The number of water molecules solvating the acid molecule were chosen such that at least three hydration shells were complete and typically ranged from 52 to 60 water molecules. All water molecules are considered explicitly in the CPMD simulations. The dimension of the simulation box was selected so as to ensure that after subtracting the van der Waals volume of the organic acid the space available for the water molecules

would

correspond to the experimental density of unity. We observed that maintaining water density as close to unity is an important consideration for obtaining accurate pKa values. The number of water molecules and the dimension of the simulation cell are different for each organic acid and are

specified in the Supporting Information, Table S1. The gradient corrected exchange

correlational DFT-HCTH/12030 functional was used in the present study because it provides significantly improved energies as compared to DFT-BLYP (Becke exchange and Lee-Yang-


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

Perr correlation) for water –ion systems.31 Dispersive interactions were accounted for by the empirical dispersion correction of Grimme.32 Core electrons were treated using the normconserving atomic pseudopotentials of Troullier-Martins,33 while valence electrons were represented in a plane wave basis set truncated at an extended energy of 70 Ry. The fictitious electron mass parameter was equal to 800 a.u.. The time step applied during the simulations was set to 6 a.u. (0.146 fs). All simulations were performed at a temperature of 300 K. The NoseHoover´ thermostat was used to control the temperature of the NVT ensemble. The initial part of the simulations was taken as the equilibration time (40 ps) and was not considered during the data analysis. The metadynamics formalism in its extended Lagrangian version was employed to explore the free-energy profile of the dissociation process as defined by the collective variable nOH, the O-H bond-distance dependent number of hydrogens coordinating the hydroxyl oxygen of the carboxylic group. =

∑

(1)

The collective variable nOH has a value close to unity in the un-dissociated state and is zero on dissociation. The instantaneous O-H bond distance is rOH and rc (= 1.6 Å) is the critical, or “cut-off”, distance beyond which the O-H bond breaks. The exponents controls the long-range behaviour and stiffness and ensures a smooth decay of the coordination number.6 The free energy profile was sampled by placing history-dependent Gaussian bias potentials along the dynamical time evolution path. The height and width of the Gaussian bias potentials were 0.0005 and 0.1 a.u., respectively, with the deposition time step 0.04 ps. The sampling was continued till


8

Page 9 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the entire free-energy profile was explored and the motion in the collective variable space was diffusive. RESULTS AND DISCUSSIONS The free-energy profile for the dissociation of acetic acid is shown in Figure 1a. These calculations involved placing the acid molecule in a cubic box of side 12.26Å along with 60 water molecules. During the metadynamics simulations dissociation is observed within 5.2 ps (see inset of Figure 1a). Two well-defined minima are observed; the deeper minimum on the left corresponds to the un-dissociated neutral molecule and that on the right to the dissociated acetate anion. In addition, a shallow minima is seen separating the reactant and product states. Snapshots along the trajectory at the locations indicated on Figure 1a are shown in Figures 1b (a movie file of configurations along the trajectory is provided in the Supporting Information, Movie S1). It may be seen from the snapshots that the first minimum (labeled 1 in Figure 1a) corresponds to the solvated neutral acetic acid. After 3 ps the proton is abstracted by a water molecule but the hydronium ion and the acetate ion remain as a bound pair; the shallow minima (labelled 2) in the free energy profile corresponds to a contact-ion pair (the geometry of the contact-ion par is described in the Supporting Information, Figure S1). Subsequently the contact-ion pair breaks-up with the hydronium ion moving away to complete the dissociation process (the second minima, labeled 3, in Figure 1a). The free-energy profiles are identical irrespective of whether the simulations were initiated from the neutral or dissociated species. It is possible to estimate the pKa values from the free energy profiles in Figure 1a as the difference in the free energy of the dissociated and un-dissociated species, pKa = ∆G / 2.303 RT. The value of pKa for acetic acid estimated from the free energy profile (Figure 1a) is 4.8; the


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 25

Figure 1: (a) The free-energy profile as a function of the collective variable, nOH, the distancedependent coordination number of the protons bound to the carboxylic group, for the dissociation of acetic acid. The difference in free energies of the dissociated and neutral species, as well as the estimated pKa value are indicated. The value in parentheses is the experimentally determined pKa value of acetic acid. The inset shows the evolution of the collective variable, nOH, during the simulation (b) Snapshots along the trajectory at the locations indicated on Figure 1a. The proton being abstracted is highlighted in yellow and also the water molecule involved at the instant of dissociation.


10

Page 11 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


value is almost identical to the experimentally determined value, 4.76. We have also examined the influence of the density and the number of water molecules on the estimates of the pKa value. The dimensions of the box are critical; the density of water should be as close to unity as possible. The densities were estimated after subtracting the van der Waals volume of acetic acid from the box volume and dividing the remainder by the product of the number of water molecules and their van der Waal volume. When the side of the cubic box is increased from 12.26 Å to 12.4Å corresponding to a water density of 0.94 g/cm3 the estimated pKa value is 5.24 and when the dimension is reduced to 12 Å corresponding to a water density of 1.04 g/cm3 the estimated pKa value is 5.04 (see Supporting Information, Figure S2). It thus appears that getting the water density as close to unity is critical to getting accurate estimates of pKa values. Although we have not examined the origin of the effect in detail it appears that lower densities makes dissociation more difficult ( at a density of 0.94 g/cm3 dissociation is observed after 7.5 ps as compared to 5.2 ps when density is unity) while higher densities favor the re-association reaction. The total number of water molecules does not appear to be as critical, provided the density is maintained as close to unity as possible. For the dissociation of acetic acid simulations using 43 and 60 water molecules returned similar estimates of the pKa value and it was only when the number was reduced to 30 water molecules that the estimate showed significant departure from the experimental value (see Supporting Information, Figure S3). Considering the nature of approximations used in this computation, it would be tempting to dismiss these estimates of the pKa value as “too good to be true”. A number of factors are responsible for the lack of accuracy in estimates of the absolute values of free energy obtained from a CPMD–metadynamics simulation. The main sources of error are the choice of the functional and secondly the errors associated with the metadynamics sampling technique.


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 25

Quantum effects as well as the zero-point energy are also ignored in these calculations.34 It may, however, be pointed out that here the pKa values are estimated from the difference in free energy values associated with the neutral and dissociated species and consequently,

there is a

cancellation of errors. Errors associated with the use and choice of functionals are also considerably reduced when one considers differences. For example, in the case of water clusters it has been reported that difference in binding energies of conformers estimated by density functionals are much closer to values obtained using ωB97X-D, LC-ωPBE-D3 (with damping to zero), and M05-2X methods as compared to absolute values of the binding energies of the conformers.35 The errors in the free-energy values introduced by the metadynamics sampling depend essentially on the sampling parameters (see previous section) and are estimated to be of the order of 1.9 kBT , for the acetic acid dissociation.36 This error would be considerably reduced when the difference in free energy values are considered. It is therefore not surprising that the estimate of the pKa value of acetic acid, obtained as the difference in free-energy values, is in good agreement with experiment. The utility of the CPMD-metadynamics formalism in predicting the dissociation constant of weak acids in general, would, however, depend on the cancellation of errors being universal, and not co-incidental to the acetic acid dissociation reaction. We show in the following sections that the method consistently predicts pKa values in good agreement with experiment for a diverse range of organic acids. One of the key features of the procedure adopted here is that all water molecules are treated explicitly at the same level of theory. Inductive Effect One of the well established factors that influences the acidity of organic acids is the inductive effect. The presence of electron withdrawing or donating groups can either stabilize or destabilize the anion and hence the magnitude of the pKa value of the acid. The effect falls off


12

Page 13 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Estimated and experimental pKa values of organic acids at 300 K, highlighting the role of the inductive effect.

pKa Acid

Calculated

Experimental *

H-COOH

3.86

3.78

CH3-COOH

4.80

4.76

CH3-CH2-COOH

4.91

4.87

CH3-CH2-CH2-COOH

4.89

4.82

ClCH2-COOH

2.92

2.86

Cl2CH-COOH

1.35

1.29

Cl3C-COOH

0.75

0.65

CH3-CH2-ClCH-COOH

2.89

2.84

CH3-ClCH-CH2-COOH

4.11

4.06

ClCH2-CH2-CH2-COOH

4.58

4.52

*from reference 37,38

rapidly with distance. We have investigated the dissociation of a series of R-COOH ( R = H,CH3, CH2CH3, CH2CH2CH3 ) acids by the CPMD-metadynamics procedure. In each case the acid molecule was placed along with water molecules in a cubic box of appropriate dimensions (see Supporting Information, Table S1) and the free-energy profile as defined by the collective variable, nOH, computed. The free-energy profiles exhibited two well-defined minima that corresponded to the neutral and dissociated acid species (see Supporting Information, Figure S4).


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 25

For all profiles a shallow minima is seen separating the reactant and product states corresponding to the formation of a contact-ion pair. The dissociation constants of the acids were estimated from the difference in the free-energy values of the dissociated and neutral species. The estimated and experimental pKa values are shown in Table 1. It may be seen that the pKa values predicted from the CPMD-metadynamics simulations are in good agreement with the experimental values. The observed trend in the pKa values with increasing R in R-COOH is clearly reproduced. The effect of electron withdrawing groups on pKa values was also investigated by considering a series of chlorine substituted aliphatic carboxylic acids. The free-energy profiles of the dissociation process as defined by the collective variable nOH were computed by the CPMDmetadynamics procedure. For all chloro-substituted acids the profiles showed two free energy minima that correspond to the neutral and dissociated species (see Supporting Information, Figure S5) and the pKa values estimated from the difference in their free energy values . The estimated values of pKa along with the experimentally determined values are shown in Table 1. It may be seen that the values are comparable. The calculated values clearly reproduce the expected influence of the number of halogen substituents as well as the effect of the distance of the electron withdrawing chlorine atom from the carboxylic acid group on acidity. Hydrogen Bonding One of the factors that strongly influences pKa values is hydrogen bonding. In acids where the dissociated anion can be stabilized by intra-molecular hydrogen bonding, the acidity is strongly enhanced. This is best illustrated by comparing pKa values of isomers, or conformers, of substituted carboxylic acids. In isomers where a hydroxyl or carboxylic group is in close proximity to the ionized carboxylic group and there is the possibility of the dissociated species being stabilized, the acidity is strongly enhanced as compared to the


14

Page 15 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: (a) The free energy profiles as function of the collective variable, nOH, for the dissociation of maleic and fumaric acids (b) Snapshots along the trajectories located close to the free energy

minima associated with the ionic dissociated states of maleic and fumaric acid. The abstracted proton is highlighted in yellow.


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 25

Table 2: Estimated and experimental pKa values of organic acids, highlighting the role of hydrogen bonding in stabilizing the dissociated acid.

pKa Calculated

Experimental *

Maleic acid

1.99

1.92

Fumaric acid

3.10

3.02

2-hydroxybenzoic acid

3.17

2.98


4.16

4.08


4.66

4.58

* from reference 37,38

isomers where the groups are spatially separated and stabilization by hydrogen bonding absent. We have investigated the dissociation of maleic (cis-butenedioic) and fumaric (transbutenedioic) by the CPMD-metadynamics procedure. The free energy profiles as defined by the collective variable nOH, for the dissociation of the two acids are shown in Figure 2a. As in the previous systems, two well defined minima associated with the neutral and ionized species are observed. The pKa values were estimated from the difference in the free energy values at the minima and are shown in Table 2 along with the experimentally determined values. The simulations correctly predict the large difference in pKa values of maleic and fumaric acids and the magnitudes, too, are comparable to the experimental values. Although the focus here was not to describe the dissociation process it was clear from a preliminary examination of the trajectories that in the case of maleic acid the ionized species is stabilized by hydrogen bonding.


16

Page 17 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3: The free energy profiles for acid dissociation of ortho-, meta- and para- isomers of hydroxybenzoic acid as a function of the collective variable nOH.

This is evident from from snapshots along the trajectories located close to the free energy minima associated with the dissociated states for maleic and fumaric acid ( Figure 2b). The dissociation of the ortho-, meta- and para- isomers of hydroxybenzoic acid were also investigated by the CPMD-metadynamics procedure. The experimental pKa values of the three isomers are very different and has been attributed to the possibility of hydrogen bonding stabilization of the ortho- hydroxybenzoate anion and its absence in the other two isomers. The computed free-energy profiles of the dissociation reaction of the hydroxyl benzoic acids are


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 25

shown in Figure 3 and the pKa values estimated from the difference in free-energy values given in Table 2. It may be seen that the simulations faithfully reproduce the experimentally observed trend and the magnitudes, too, are comparable to the experimental values.

CONCLUSIONS In conclusion we have shown that quantum CPMD based metadynamics simulations can predict pKa values of weak acids with reasonable accuracy. The formalism treats all solvent water molecules explicitly while retaining conceptual simplicity. The number of water molecules considered was sufficient to complete the three hydration shells surrounding the acid molecule. We use the distance-dependent coordination number of the protons bound to the hydroxyl oxygen of the carboxylic group as the collective variable to explore the free energy profile of the dissociation process. Two distinct minima corresponding to the dissociated and un-dissociated states of the acid are observed and the difference in their free energy values provides the estimate for pKa, the acid dissociation constant. Since the pKa values are estimated from the difference in free-energy values the procedure benefits from a reduction in errors due to cancellation, which appears to hold true for all systems studied and suggests that this may be universally valid. We show that the method predicts the pKa value of acetic acid in good agreement with experiment and then go on to show for a series of organic acids that ab initio CPMD simulations in conjunction with metadynamics calculations of the free energy profile of the dissociation reaction can provide reasonable estimates of the pKa value. The acids investigated


18

Page 19 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


were aliphatic carboxylic acids, chlorine substituted carboxylic acids, cis- and transbutenedioic, and the isomers of hydroxy-benzoic acid. These systems were chosen to highlight whether the procedure could correctly account for the influence of the inductive effect as well as hydrogen bonding on pKa values of weak organic acids. In both situations the CPMD-metadynamics procedure faithfully reproduces the experimentally observed trend and the magnitudes of the pKa values, too, are comparable to the experimentally determined values. The results presented here clearly validate the usefulness of this procedure in estimating pKa values and as we show, in subsequent studies, is easily extended to larger systems like amino acids and peptides. ASSOCIATED CONTENT Supporting Information Available: (Table S1) The number of water molecules and dimension of the simulation cell for the acid molecules investigated. (Movie S1) Movie file of the configurations along the trajectory of the dissociation of acetic acid. (Figure S1) Structure of the contact-ion pair formed during the dissociation of Acetic acid. (Figure S2) The effect of water density on pKa values. (Figure S3) The effect of the number of water molecules on pKa values. (Figure S4) The free energy profiles for the dissociation of a series of R-COOH (R = H,CH3, CH2CH3, CH2CH2CH3) (Figure S5) The free energy profiles for the dissociation of a series of chlorine substituted carboxylic acids. This material is available free of charge via the Internet at http://pubs.acs.org.

ACKNOWLEDGMENTS The authors thank the Supercomputer Education and Research Centre at the Indian Institute of Science, Bangalore for computational facilities.


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 25

REFERENCES (1) Alongi, K. S.; Shields, G. C. Theoretical Calculations of Acid Dissociation Constants: A Review Article. Annu. Rep. Comput. Chem. 2010, 6, 113-138. (2) Alexov, E.; Mehler, E. L.; Baker, N.; Baptista, A. M.; Huang, Y.; Milletti, F.; Nielsen, J. E.; Farrell, D.; Carstensen, T.; Olsson, M. H.; Shen, J. K.; Warwicker, J.; Williams, S.; Word, J. M. Progress in the Prediction of pKa Values in Proteins. Proteins 2011, 79, 3260-3275. (3) Ho, J.; Coote, M. L. A Universal Approach for Continuum Solvent pKa Calculations: Are we there Yet? Theor. Chem. Acc. 2010, 125, 3-21. (4) Car, R.; Parrinello, M.

Unified Approach for Molecular Dynamics and Density-

Functional Theory. Phys. Rev. Lett. 1985, 55, 2471-2474. (5) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. PNAS 2002, 99, 12562-12566. (6) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302-4. (7) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 826. (8) Hassanali, A. A.; Cuny, J.; Verdolino, V.; Parrinello, M. Aqueous Solutions: State of the Art in Ab initio Molecular Dynamics. Philos. Trans. R. Soc. A 2014, 372, 20120482. (9) Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. Absolute pKa Determinations for Substituted Phenols. J. Am. Chem. Soc. 2002, 124, 6421-6427.


20

Page 21 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10) Toomsalu, E.; Koppel, I. A.; Burk, P. Critical Test of Some Computational Chemistry Methods for Prediction of Gas-Phase Acidities and Basicities. J. Chem. Theory Comput. 2013, 9, 3947−3958. (11) Lim, C.; Bashford, D.; Karplus, M. Absolute pKa Calculations with Continuum Dielectric Methods. J. Phys. Chem. 1991, 95, 5610-5620. (12) 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. (13) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion−Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066-16081. (14) 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, 11331152. (15) Pliego, J. R.; Riveros, J. M. The Cluster-Continuum Model for the Calculation of the Solvation Free Energy of Ionic Species. J. Phys. Chem. A 2001, 105, 7241−7247. (16) Pliego, J. R., Jr; Riveros, J. M. Theoretical Calculation of pKa Using the ClusterContinuum Model. J. Phys. Chem. A 2002, 106, 7434-7439.


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 25

(17) 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. (18) Ho, J. Predicting pKa in Implicit Solvents: Current Status and Future Directions. Aust. J. Chem. 2014, http://dx.doi.org/10.1071/CH14040. (19) Choi, C. H.; Re, S.; Feig, M.; Sugita, Y. Quantum Mechanical/Effective Fragment Potential Molecular Dynamics (QM/EFP-MD) Study on Intra-Molecular Proton Transfer of Glycine in Water. Chem. Phys. Lett. 2012, 539, 218−221. (20) Ghosh, M. K.; Re, S.; Feig, M.; Sugita, Y.; Choi, C. H. Interionic Hydration Structures of NaCl in Aqueous Solution: A Combined Study of Quantum Mechanical Cluster Calculations and QM/EFPMD Simulations. J. Phys. Chem. B 2013, 117, 289−295. (21) Ghosh, M. K.; Uddin, N.; Choi, C. H. Hydrophobic and Hydrophilic Associations of a Methanol Pair in Aqueous Solution. J. Phys. Chem. B 2012, 116, 14254−14260. (22) Uddin, N.; Choi, T. H.; Choi, C. H. Direct Absolute pKa Predictions and Proton Transfer Mechanisms of Small Molecules in Aqueous Solution by QM/MM-MD. J. Phys. Chem. B 2013, 117, 6269−6275. (23) Nelson, J. G.; Peng, Y.; Silverstein, D. W.; Swanson, J. M. J. Multiscale Reactive Molecular Dynamics for Absolute pKa Predictions and Amino Acid deprotonation. J. Chem. Theory Comput. 2014, 10, 2729−2737. (24) Park, J.M.; Laio, A.; Iannuzzi, M.; Parrinello, M. Dissociation Mechanism of Acetic Acid in Water. J. Am. Chem. Soc. 2006, 128, 11318-11319.


22

Page 23 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(25) Galib, M.; Hanna, G. Mechanistic Insights into the Dissociation and Decomposition of Carbonic Acid in Water via the Hydroxide Route: An Ab Initio Metadynamics Study. J. Phys. Chem. B. 2011, 115, 15024-15035. (26) Lee, J. G.; Asciutto, E.; Babin, V.; Sagui, C.; Darden, T.; Roland, C. Deprotonation of Solvated Formic Acid: Car-Parrinello and Metadynamics Simulations. J. Phys. Chem. B. 2006, 110, 2325-2331. (27) Ivanov, I.; Chen, B.; Raugei, S.; Klein, M. L. Relative pKa Values from First-Principles Molecular Dynamics: The Case of Histidine Deprotonation. J. Phys. Chem. B 2006, 110, 6365−6371. (28) Sulpizi, M.; Sprik, M. Acidity Constants from Vertical Energy Gaps: Density Functional Theory Based Molecular Dynamics Implementation. Phys. Chem. Chem. Phys. 2008, 10, 5238–5249. (29) CPMD, version 3.13.2; IBM Corp 1990−2008, MPI für Festkörper-forschung Stuttgart 1997−2001. http://www.cpmd.org/.

(30) Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M. New Generalized Gradient Approximation Functionals. J. Chem. Phys. 2000, 112, 1670-1678. (31) Izvekov, S.; Voth, G. A. Ab Initio Molecular Dynamics Simulation of Aqueous Proton Solvation and Transport Revisited. J. Chem. Phys. 2005, 123, 044505-9. (32) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a LongRange Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799.


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 25

(33) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for Plane-Wave Calculations. Phys. Rev. B 1991, 43, 1993-2006. (34) Garcia-Viloca, M.; Alhambra, C.; Truhlar, D. G.; Gao, J. Inclusion of QuantumMechanical Vibrational Energy in Reactive Potentials of Mean Force. J. Chem. Phys. 2001, 114, 9953-9958. (35) Leverentz, H. R.; Qi, H. W.; Truhlar, D. G. Assessing the Accuracy of Density Functional and Semiempirical Wave Function Methods for Water Nanoparticles: Comparing Binding and Relative Energies of (H2O)16 and (H2O)17 to CCSD(T) Results. J. Chem. Theory Comput. 2013, 9, 995−1006. (36) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics. J. Phys. Chem. B 2005, 109, 6714-6721. (37) Braude, E. A.; Nachod, F. C. Determination of Organic Structures by Physical Methods. New York: Academic Press, 1955. (38) Perrin, D.D. Dissociation Constants of Organic Bases in Aqueous Solution. Butterworths, 1972.


24

Page 25 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents


25

Estimating successive pKa values of polyprotic acids from ab initio molecular dynamics using metadynamics: the dissociation of phthalic acid and its isomers.

N2 dissociation on W(110): An ab initio molecular dynamics study on the effect of phonons.

The hydrogen diffusion in liquid aluminum alloys from ab initio molecular dynamics.

Reactive molecular dynamics models from ab initio molecular dynamics data using relative entropy minimization.

Ab initio metadynamics study on hydronium ion dynamics at acid-functionalized interfaces: effect of surface group density.

Li adsorption, hydrogen storage and dissociation using monolayer MoS2: an ab initio random structure searching approach.

Ab initio molecular dynamics study of the aqueous HOO(-) ion.

Resolving the anomalous infrared spectrum of the MeCN-HCl molecular cluster using ab initio molecular dynamics.

Dissociation and recombination of D₂ on Cu(111): ab initio molecular dynamics calculations and improved analysis of desorption experiments.

mineral interfaces from ab initio molecular dynamics simulations.

Dissolution of NaCl nanocrystals: an ab initio molecular dynamics study.

Ab initio molecular dynamics study of small alkali metal clusters.

Simulations of the dissociation of small helium clusters with ab initio molecular dynamics in electronically excited states.

Predicting pKa Values in Aqueous Solution for the Guanidine Functional Group from Gas Phase Ab Initio Bond Lengths.

Multiple time step integrators in ab initio molecular dynamics.

Structural and electronic properties of organo-halide hybrid perovskites from ab initio molecular dynamics.

Halogen and Hydrogen Bonding Benzothiophene Diol Derivatives: A Study Using ab initio Calculations and X-Ray Crystal Structure Measurements.

Benchmarking ab initio binding energies of hydrogen-bonded molecular clusters based on FTIR spectroscopy.

Slide fastener reduction of graphene-oxide edges by calcium: insight from ab initio molecular dynamics.

Liquid Aluminum: atomic diffusion and viscosity from ab initio molecular dynamics.

Raman spectra from ab initio molecular dynamics and its application to liquid S-methyloxirane.

Ab initio and metadynamics studies on the role of essential functional groups in biomineralization of calcium carbonate and environmental situations.

Ab initio molecular dynamics of solvation effects on reactivity at electrified interfaces.

Investigation of the hydroxylation mechanism of noncoupled copper oxygenases by ab initio molecular dynamics simulations.