Microscopic description of acid–base equilibrium

Grifoni, Emanuele; Piccini, GiovanniMaria; Parrinello, Michele

doi:10.1073/pnas.1819771116

Microscopic description of acid–base equilibrium

Emanuele Grifoni, GiovanniMaria Piccini, and Michele Parrinello [email protected]Authors Info & Affiliations

Contributed by Michele Parrinello, January 10, 2019 (sent for review November 29, 2018; reviewed by Christoph Dellago and David E. Manolopoulos)

February 14, 2019

116 (10) 4054-4057

https://doi.org/10.1073/pnas.1819771116

PDF/EPUB

Significance

Acid–base reactions are among the most important chemical processes. Yet we lack a simple way of describing this class of reactions as a function of the atomic coordinates. In fact, once dissolved in water,

H^{+}

and its conjugate anion

{OH}^{-}

have a highly fluxional structure difficult to pin down. Here we solve this issue by taking the point of view of describing acid–base reactions as an equilibrium between the solute and the whole solvent. This allows identifying generally applicable descriptors. As a consequence it is now possible to perform quantitative enhanced sampling simulation of acid–base reaction in water and in other environments such as the zeolite cavities or at surfaces.

Abstract

Acid–base reactions are ubiquitous in nature. Understanding their mechanisms is crucial in many fields, from biochemistry to industrial catalysis. Unfortunately, experiments give only limited information without much insight into the molecular behavior. Atomistic simulations could complement experiments and shed precious light on microscopic mechanisms. The large free-energy barriers connected to proton dissociation, however, make the use of enhanced sampling methods mandatory. Here we perform an ab initio molecular dynamics (MD) simulation and enhance sampling with the help of metadynamics. This has been made possible by the introduction of descriptors or collective variables (CVs) that are based on a conceptually different outlook on acid–base equilibria. We test successfully our approach on three different aqueous solutions of acetic acid, ammonia, and bicarbonate. These are representative of acid, basic, and amphoteric behavior.

Acid–base reactions play a key role in many branches of chemistry. Inorganic complexation reactions, protein folding, enzymatic processes, polymerization, catalytic reactions, and many other transformations in different areas are sensitive to changes in pH. Understanding the pH role in these reactions implies having control over their reactivity and kinetics.

The crucial importance of pH has stimulated the collection of a large amount of data on acid–base equilibria. These are typically measured in gas and condensed phases, using spectroscopic and potentiometric techniques. However, there are practical limitations to the accuracy of these methods especially in condensed phases (1). Furthermore it is very difficult to extract from experimental data a microscopic picture of the processes involved. It is thus not surprising that acid–base equilibrium has been the subject of intense theoretical activity (1–12).

The acidity of a chemical species in water can be expressed in terms of

p K_{a}

, the negative logarithm of the acid dissociation constant. There are two ways of calculating these values, one static and the other dynamic.

The most standard approach is the static one in which solution-phase free energies, and consequently

p K_{a}

s, are obtained by closing a Born–Haber cycle composed of gas phase and solvation free energies (1, 3–7). While extremely successful in many cases, the static approach has some limitations. A solvation model needs to be chosen and continuum solvent models have a limited accuracy. This is particularly true in systems like zeolites or proteins characterized by irregular cavities in which an implicit description of the solvent is challenging. Obviously from such an approach dynamic information cannot be gained. Furthermore, there can be competitive reactions that cannot be taken into account unless explicitly included in the model.

In principle these limitations could be lifted in a more dynamical approach based on molecular dynamics (MD) simulations in which the solvent molecules are treated explicitly. If one had unlimited computer time, such simulations would explore all possible pathways and assign the relative statistical weight to the different states. Unfortunately the presence of kinetic bottlenecks frustrates this possibility by trapping the system in metastable states, since different protonation states are separated by large barriers. Furthermore in acid–base reactions chemical bonds are broken and formed. This requires the use of ab initio MD in which the interatomic forces are computed on the fly from electronic structure theories. This makes the calculation more expensive and reduces further the time scale that can be explored.

To overcome this difficulty, the use of enhanced sampling methods (13) that accelerate configurational space exploration becomes mandatory. A very popular class of enhanced sampling methods is based on the identification of the degrees of freedom that are involved in the slow reaction of interest. These degrees of freedom are usually referred to as collective variables (CVs) and are expressed as explicit functions of the atomic coordinates

R

. Sampling is then enhanced by adding a bias that is a function of the chosen CVs (14–16). Furthermore, designing a proper set of good CVs has also a deeper meaning. Successful CVs capture in a condensed way the physics of the problem, identify its slow degrees of freedom, and lead to a useful modelistic description of the process.

In standard chemical reactions, this is relatively simple since well-defined structures can be assigned to reactants and products (17–19). This is not the case for acid–base reactions in which a proton is added to or subtracted from the solute. Once this process has taken place, water ions (

H^{+}

or/and

{OH}^{-}

) are solvated and their structure becomes elusive. In fact, water ions can rapidly diffuse in the medium via a Grotthuss mechanism (20). They became highly fluxional and the identity of the atoms taking part in their structure changes continuously. The nature of these species is thus difficult to capture in an explicit analytic function of

R

. However, given the relevance of acid–base reactions, many attempts have been made at defining these entities (8–12). Unfortunately these CVs have an ad hoc nature and, while successful in this or that case, cannot be generally applied.

To build general and useful CVs we make two conceptual steps. One is to look at the acid–base process as a reaction involving only a few moieties, namely the whole solvent and the reacting residues in the solvated molecule. For example, when there is only one type of dissociating residue, we think of the acid–base equilibrium as a reaction of the type

A + H_{2 N} O_{N} ⇌ B^{q_{0}} + H_{2 N + q_{1}} O_{N}^{q_{1}},

[1]

where

N

is the number of water molecules,

A

and

B

are a generic acid–base molecule in solution and its conjugate species, respectively,

q_{0}

and

q_{1}

are integers that can assume values

+ 1

and −1 according to the acid–base behavior of the species, and

q_{1} + q_{0} = 0

.

This implies that we do not look at the solvent as a set of molecules that compete to react with the acid–base species. Rather we consider the solvent in its entirety as one of the two adducts. Taking this point of view is especially relevant in polar solvents like water that are characterized by highly structured networks. In this case the presence of an excess or a deficiency of protons changes locally the network structure and this distortion propagates along the entire network.

Since the very early days of Wicke and Eigen (21) and Zundel and Metzger (22), researchers have struggled with how many molecules should be included in the definition of the perturbation (23–25). Given the absence of physical parameters capable of giving a clear and unequivocal answer to this question, the idea of considering the solvent as a whole circumvents this problem. Thus, the solvent is not just a medium with a passive role, but it is looked at as an ensemble of molecules that contribute collectively to the formation of the conjugate acid–base pair. This point of view is much closer to the original one proposed by Brønsted and Lowry in which the reaction can be seen as a simple exchange of a hydrogen cation between an acid–base pair.

For the reaction to take place the center of the perturbation has to move away from the solute. Thus, the second important step is to monitor the center of the perturbation. Due to Grotthuss-like mechanisms, the perturbation moves along the network. This can lead to different definitions of the defect center. However, if we tessellate the whole space using Voronoi polyhedra centered on water oxygen atoms, we can assign unequivocally every hydrogen atom to one and only one of these polyhedra. The site whose Voronoi polyhedron contains an anomalous number of protons is taken as the center of the perturbation (Fig. 1).

Fig. 1.

Two examples of partitioning the space. (*Left*) We show a convectional approach in which the distance from the oxygen atom is used to define its surroundings. Clearly artificial superpositions can be seen. (*Right*) The Voronoi tessellation does not suffer from these shortcomings.

This point of view gives the method a very general nature, making it applicable to every acid–base system, without the need of fixing beforehand the reacting pairs. Thus, it is possible to explore all of the relevant protonation states even in systems composed of more than one acid–base pair.

This general approach allows defining CVs without having to impose specific structures or select the identity of the atoms involved. We test our method by performing metadynamics simulations in a weak acid case (acetic acid), in a weak base (ammonia), and in an amphoteric species (bicarbonate) chosen as benchmarks because of their comparable strength, but different acid–base behavior.

Methods

As discussed above we introduce two CVs, one related to the protonation state and another that locates the charge defects and measures their relative distance. Both of these CVs need a robust definition for assigning the hydrogen atoms to the respective acid–base site. To achieve this result we partition the whole space into Voronoi polyhedra centered on the acid–base sites

i

located at

R_{i}

. The sites include all of the atoms able to break and form bonds with an acid proton. The standard Voronoi space partition is described by a set of index functions

w_{i} (r)

centered on the different

R_{i}

s such that

w_{i} (r) = 1

if the

i

th atom is the closest to

r

and w_i(r) = 0 otherwise. For their use in enhanced sampling methods CVs need to be differentiable. To this effect we introduce a smooth version of the index functions,

w_{i}^{s} (r)

. These are defined using softmax functions

w_{i}^{s} (r) = \frac{e^{- λ | R_{i} - r |}}{\sum_{m} e^{- λ | R_{m} - r |}},

[2]

where

i

and

m

run over all the acid–base sites and

λ

controls the steepness with which the curves decay to 0, that is, the selectivity of the function. With an appropriate choice of

λ

this definition achieves the desired result as shown in Fig. 2. In such a way, a hydrogen atom in a position

R_{j}

is assigned to the polyhedron centered on the site

i

with the weight

w_{i} (R_{j})

. Then, the total number of hydrogen atoms assigned to the

i

th acid–base site is

W_{i} = \sum_{j \in H} w_{i}^{s} (R_{j}),

[3]

where the summation on

j

runs over all the hydrogen atoms.

Fig. 2.

Smooth tessellation of a 2D space with cells centered on the three water molecule oxygen atoms. The flat blue regions represent the portion of space in which the function assumes a value of 1 and the yellow ones represent the borders among cells. This surface has been obtained with a value of $λ = 4$ .

One can associate to each acid–base site a reference value

W_{i}^{0}

that counts the number of bonded hydrogen atoms in the neutral state. The difference between the instantaneous value of hydrogen atoms and the reference one is

δ_{i} = W_{i} - W_{i}^{0} .

[4]

When different from zero,

δ_{i}

will signal whether the

i

th site has gained or lost a proton. In the case of water oxygen atoms, a hydronium ion has a

δ_{i} = + 1

while a hydroxide ion has

δ_{i} = - 1

.

We then group the acid–base sites in species. For instance, in the case of the simplest amino acid glycine in aqueous solution the number of species

N^{s}

will be equal to 3. All water oxygen atoms belong to one species; then one counts in another species the two carboxylic oxygen atoms, and finally one considers as the third species the nitrogen atom of the amino group.

In the spirit of this work we count the total excess or defects of protons associated to each species:

q_{k} = \sum_{i \in k} δ_{i} .

[5]

This implies that we are not interested in the specific identity of the reacted site, but whether or not the

k

th species in its entirety has increased (

q_{k} = + 1

), has decreased (

q_{k} = - 1

), or has not changed its number of protons. If we consider a solute with only one reactive moiety, then each possible state of the system can be described by one of the three 2D vectors (0,0), (−1,1), or (1,−1).

In the general case each protonation state can be described by a vector

\vec{q} = (q_{0}, q_{1}, \dots q_{N^{s} - 1})

with dimension equal to the number of inequivalent reactive sites,

N^{s}

. A more exhaustive explanation is provided in SI Appendix.

For use in enhanced sampling these vectors need to be expressed as a scalar function

f = f (\vec{q})

such that, for each physically relevant

\vec{q}

,

f

attains values able to distinguish the different overall protonation states. There are infinitely many ways of constructing a scalar from a vector. Possibly the simplest choice is to write

f (\vec{q}) = \vec{X} \cdot \vec{q}

and, to distinguish between different protonation states, to choose

\vec{X} = (2^{0}, 2^{1}, 2^{2}, \dots 2^{N^{s} - 1})

.

This leads to the following definition for the CV that is used to describe the protonation state of the system,

s_{p} = \sum_{k = 0}^{N^{s} - 1} 2^{k} \cdot q_{k},

[6]

where

k

are the indexes used to label the respective reactive site groups. In SI Appendix an example is worked out in detail. Of course the CV is made continuous by the use of

W_{i}

in the calculation of the

δ_{i}

needed to evaluate

q_{k}

in Eq. 5.

The second CV is a summation of distances between all acid–base sites multiplied by their partial charge

δ_{i}

,

s_{d} = \sum_{i, m > i} - r_{i m} \cdot δ_{i} \cdot δ_{m},

[7]

where the indexes

i

and

m

run over all the acid–base sites belonging to different

k

groups, and

r_{i m}

is the distance between the two atoms. In this way, just the acid–base pair that has exchanged a proton gives a contribution different from zero. Eq. 7 is valid only when one single-conjugate acid–base pair is present. However, due to the action of bias, it may occur occasionally that several acid–base pairs are formed. To avoid sampling these very unlikely events we apply a restraint on the number of pairs. Further details are provided in SI Appendix.

Results

We have applied our method to three aqueous solutions of acetic acid, ammonia, and bicarbonate as representations of a weak acid, a weak base, and an amphoteric compound, respectively. The setups of all three simulations are identical except for the identity of the solvated molecules. This ensures that the outcome reflects the different chemistry of these three systems and that there is no bias due to the initial condition.

Each simulation of the systems was performed with Born–Oppenheimer MD simulations combined with well-tempered metadynamics (14, 26) using the CP2K package (27) patched with PLUMED 2 (28) and strongly constrained and appropriately normed functional (29) for the xc energy,

E_{x c}

. See SI Appendix for details.

In Fig. 3 we plot the free-energy surfaces (FESs) as a function of

s_{p}

and

s_{d}

. These FESs vividly reproduce the expected behavior. They all have a minimum at

s_{p} = 0

that corresponds to the state in which no charges are present in the solvent. In the acetic acid FES (Fig. 3A) a second minimum close to

s_{p} = - 1

reflects its acid behavior. By contrast, the ammonia FES (Fig. 3B) shows a second minimum close to

s_{p} = 1

. The shapes of ammonia and acetic acid FESs are approximately related by a mirror symmetry reflecting their contrasting behavior. Similarly the bicarbonate symmetric FES (Fig. 3C) mirrors its amphoteric character.

Fig. 3.

(*A–C*) Free-energy surfaces along $s_{p}$ and $s_{d}$ of acetic acid (A), ammonia (B), and bicarbonate (C) in aqueous solution. Color bars indicate the free energy expressed in kJ $\cdot$ mol⁻¹ units. The CV $s_{d}$ is expressed in angstroms.

As the conjugate pair is formed

s_{d}

starts to assume positive values corresponding to the separation and diffusion of the conjugate pair. Compared with the undissociated state in which only

s_{d} = 0

is allowed, states where a conjugate pair is present show an elongated shape of the basins along this variable. This is caused by the diffusive behavior of the hydronium and hydroxide ions in solution that makes accessible a continuum range of distances. Moreover, along this CV we can observe a barrier around 1.5 corresponding to the breaking of the covalent bond between the hydrogen atom and the acid–base site.

Conclusions

The general applicability of this method to systems with different natures is an important step made in their understanding and description. The scheme can be extended to include quantum nuclear effects with the use of path integral molecular dynamics (30). This would be of quantitative significance since for instance

p K_{a}

values are affected by deuteration. Moreover, the absence of assumptions or impositions about reactive candidates or reaction paths allows extending this method to systems of increasing complexity which cannot be addressed with traditional methods. Examples of questions that can now be answered are tautomeric equilibria in biochemical processes and acid behavior in zeolites and on the surface of oxides exposed to water.

Acknowledgments

Calculations were carried out on the ETH Euler cluster and on the Mönch cluster at the Swiss National Supercomputing Center. This research was supported by the European Union Grant ERC-2014-AdG-670227/VARMET.

Supporting Information

Appendix (PDF)

Download
2.34 MB

References

1

J Ho, ML Coote, First-principles prediction of acidities in the gas and solution phase. Wiley Interdiscip Rev Comput Mol Sci 1, 649–660 (2011).

Featured Topics

Articles By Topic

Featured Topics

Articles By Topic

Featured Topic

Articles By Topic

Significance

Abstract

Sign up for PNAS alerts.

Methods

Results

Conclusions

Acknowledgments

Supporting Information

References

Information

Published in

Classifications

Copyright

Submission history

Keywords

Acknowledgments

Authors

Affiliations

Notes

Competing Interests

Metrics

Citation statements

Altmetrics

Citations

Cite this article

Cited by

View options

PDF format

Login options

Recommend to a librarian

Purchase options

Restore content access

Figures

Tables

Other

Share

Share article link

Share on social media

Further reading in this issue

Three-dimensional optical trapping and orientation of microparticles for coherent X-ray diffraction imaging

Fibroblast growth factor receptor influences primary cilium length through an interaction with intestinal cell kinase

E-protein regulatory network links TCR signaling to effector Treg cell differentiation

Bodily maps of emotions

Effect of a giant meteorite impact on Paleoarchean surface environments and life

Physician–patient racial concordance and newborn mortality

Sign up for thePNAS Highlights newsletter