Local initiation conditions for water autoionization
See allHide authors and affiliations
Edited by Phillip L. Geissler, University of California, Berkeley, CA, and accepted by Editorial Board Member R. D. Levine April 5, 2018 (received for review August 9, 2017)

Significance
The dissociation of water is arguably the most fundamental chemical reaction occurring in the aqueous phase. Despite that the splitting of a water molecule very seldom occurs, the reaction is of major importance in many areas of chemistry and biology. Direct experimental probing of the event is still impossible and also simulating the event via accurate computer simulations is challenging. Here, we achieved the latter via specialized rare-event algorithms estimating rates of dissociation in agreement with indirect experimental measurements. Even more interestingly, by a rigorous analysis of our results we identified anomalies in the water structure that act as initiators of the reaction, a finding that suggests paradigms for steering and catalyzing chemical reactions.
Abstract
The pH of liquid water is determined by the infrequent process in which water molecules split into short-lived hydroxide and hydronium ions. This reaction is difficult to probe experimentally and challenging to simulate. One of the open questions is whether the local water structure around a slightly stretched OH bond is actually initiating the eventual breakage of this bond or whether this event is driven by a global ordering that involves many water molecules far away from the reaction center. Here, we investigated the self-ionization of water at room temperature by rare-event ab initio molecular dynamics and obtained autoionization rates and activation energies in good agreement with experiments. Based on the analysis of thousands of molecular trajectories, we identified a couple of local order parameters and show that if a bond stretch occurs when all these parameters are around their ideal range, the chance for the first dissociation step (double-proton jump) increases from
Among all possible chemical reactions that occur in water, the most fundamental is the water dissociation reaction (1), which is of major importance in many areas of chemistry and biology (2). Water plays an important role as a universal solvent for a wide variety of chemical processes and can act both as an acid and as a base. In aqueous solution, water will self-ionize and form hydroxide (OH−) and hydronium (H3O+) ions which take on Eigen- or Zundel-like structures (2⇓⇓⇓–6). Experiments show that the mean lifetime for an individual molecule before undergoing autoionization is about 11 h (7, 8).
The autoionization event has not been directly probed by experiments and the dissociation rate is obtained using the water dissociation equilibrium constant and the rate for the much faster recombination reaction (7, 8). The experimental challenges make the autoionization event a pertinent target for computer simulations for which previous constrained ab initio simulations have given important information about the mechanism (9⇓–11). However, the use of constraints leads to a loss of the spontaneous dynamics of the system and the selection of a reaction coordinate that accurately measures the progress of the reaction is challenging. These limitations can be avoided by path-sampling methods such as transition path sampling (TPS) (12) or replica exchange transition interface sampling (RETIS) (13, 14) which are specifically designed for sampling rare events without altering the dynamics while less influenced by the choice of the order parameter (15). Geissler et al. (16) applied TPS with ab initio molecular dynamics (MD) to simulate just 10 uncorrelated autoionization events and demonstrated that the mechanism involves transfer of protons along a hydrogen bond wire with concomitant breaking of the wire. In their work, local solvent properties (e.g., ion coordination numbers and the presence of specific hydrogen bonds) were used to interpret the destabilization that leads to ionization. The absence of clear visually observable correlations led to the conclusion that the destabilization is caused by rare electric-field fluctuations which arise primarily from long-range electrostatic interactions and thus that local order parameters are not suitable to describe the event. Hassanali et al. (17) studied the reverse recombination reaction (i.e., neutralization of ionized water molecules) with standard ab initio MD and reported that this event takes place by a collective compression of the water wire bridging the ions, followed by a triple concerted proton jump. The OH− ion which is neutralized remains in a hypercoordinated state and Hassanali et al. (17) hypothesized that it could serve, together with the compression of the wire, as a nucleation site for autoionization. This view opposes the statement of Geissler et al. (16) that the dissociation event is primarily triggered by nonlocal structural fluctuations. We note that concerted proton transfers and collective compression of water wires have also been observed for the recombination of a weak base in water (18).
Both of these studies give important information about the autoionization mechanism, although they do not unambiguously reveal the conditions that need to accompany a bond stretch fluctuation to initiate the reaction. In this work, we aim to tackle this ambiguity and quantitatively identify initiation conditions for water autoionization. Simulating the dissociation events may not be sufficient as the apparent initiation conditions observed in trajectories that lead to dissociation may also be present in trajectories with an initial bond stretch but still fail to dissociate. Also nonreactive or “almost reactive” trajectories contain important information as these allow for identification of effective initiation conditions that really matter: those that discriminate between reactive and nonreactive trajectories. To collect this information, we applied the RETIS method and harvested reactive and nonreactive trajectories which we analyzed using the recently developed predictive power method (19) and we built a predictive machine-learning model (20). This allowed us to quantitatively examine the importance of local order parameters and initiation conditions for water autoionization. Based on this analysis we identify important initiation triggers and calculate the full rate of dissociation.
Results and Discussion
The autoionization event was investigated using ab initio RETIS simulations as described in Materials and Methods. For the RETIS simulations, we used a relatively simple geometric distance order parameter, λ, as illustrated in Fig. 1: When the system consists of only H2O species, λ is the largest covalent O–H bond distance, and when the system contains OH− and H3O+ species, λ is taken as the shortest distance between the oxygen in OH− and the hydrogen atoms in H3O+. In the following, we refer to the oxygen atom used for the order parameter as Oλ. The type of species (OH−, H2O, or H3O+) was identified by allocating to each hydrogen a single bond connecting it to the closest oxygen. Note that the definition of the order parameter does not require a threshold for defining a chemical bond nor does it constrain the order parameter to specific water molecules for the duration of the simulation. This means that we compute the rate of dissociation of any water molecule in the system instead of a single targeted O–H bond or water molecule.
The order parameter and the probability for water autoionization. (A and B) Definition of the order parameter (λ, dashed line), taken as the largest covalent O–H distance in the system when no ionic species are present (A) or as the shortest distance between the OH− oxygen atom and the hydrogen atoms in H3O+ when ionic species are present (B). A hydrogen bond wire with four members is shown with red (oxygen) and white (hydrogen) spheres and the distances
From our RETIS simulations, the water dissociation rate constant,
In the present case, we placed the final interface beyond the maximum distance obtainable in our system. All trajectories were thus propagated until the system contained only H2O species again. Separated ions may still recombine fast (within a few femtoseconds) even if the separation is large (16) and this observation was confirmed in our analysis (SI Appendix, Fig. S1). To better identify and distinguish the metastable ionized states, we used path reweighting (21) to project the crossing probability on an alternative order parameter,
In Fig. 1, we show the calculated crossing probability from our simulations as a function of the order parameter. In principle, there are two potential mechanisms which lead to an increase of the reaction coordinate λ after the first proton jump. The ionic species can separate further by another proton jump, the so-called Grotthuss mechanism, reassigning the H3O+ or OH− ion to another oxygen and causing a sudden discontinuous increase in the reaction coordinate. A second possible mechanism keeps the first ionic species intact and lets them move away from each other by diffusion, yielding a more gradual increase of the reaction coordinate. Based on the completely flat intermediate plateau region between 1.5 Å and 3.2 Å, we can conclude that only the first mechanism is effective. For
Comparing with experimentally determined dissociation constants at 25 °C [
We also calculated the average energy of the generated trajectories as a function of the order parameter (Fig. 1). The energy is expected to converge to the activation energy as can be derived from the temperature derivative of the rate constant (24, 25). We note that this activation energy gives a more direct comparison with experiments than free energy barriers which depend on the choice of order parameter. The activation energy obtained from the average energy of the accepted paths is approximately 17.8 kcal/mol. For comparison, an Arrhenius plot of the experimental data of Natzle and Moore (8) results in an activation energy of approximately 17.3 kcal/mol while Eigen and Maeyer (7) reported an activation energy between 15.5 kcal/mol and 16.5 kcal/mol. The deviation with our result is lower than the typical error margin mentioned above and the fact that the experimental activation barriers are lower than our simulation result, despite having lower rate constants, is rather remarkable. Since experimental data on this topic are at least three decades old, we hope that our finding will encourage future experimental investigations on the dissociation reaction.
Path sampling methods generate reactive (and nonreactive) trajectories which can be used to discover possible mechanisms and initiation conditions. To characterize these conditions, we considered additional collective variables, which we label
The first collective variable we consider is the length of the hydrogen bond wire bridging the nascent ion species. Our aim is to predict the outcome of initiated trajectories and in particular the initiation conditions for reactive events. Thus, we cannot define the hydrogen bond wires as connecting the ionic species, since this is one of the outcomes we wish to predict. For a single trajectory, we define the hydrogen bond wire as the shortest wire containing the Oλ species and
In addition, we have considered the following four collective variables which describe the local structure surrounding the Oλ species: (i) The number of hydrogen bonds accepted,
After defining the extra ξs, we analyzed the trajectories using the predictive power method (19). This method begins by classifying the trajectories as reactive or nonreactive based on two thresholds
We first investigated the lengths of hydrogen bond wires containing three, four, and five water molecules. Comparing the predictive abilities for these collective variables (respectively,
The concerted behavior of the autoionization event. (A) The distances (
Comparing the additional collective variables (SI Appendix, Figs. S6 and S7), we find that
Increasing the predictive power for water autoionization by considering additional collective variables. (A) The predictive power (
Inspecting the initiation conditions in more detail, we investigate the reactive and nonreactive distributions
Initiation conditions and local collective variables. (A) Reactive (
If we consider the q coordinate, we observe that
Machine learning (ML) applied to path-sampling data (33, 34) is a promising approach to find important collective variables that can easily be missed by human intuition. To explore this possibility, we built ML models for predicting the outcome of trajectories given the state of the water system early in the trajectories. We focus on the same range as in the predictive power analysis and we use the state of the system, when
In addition, to avoid a potential risk of overinterpretation we opted to restrict the complexity of the ML decision process and imposed a maximum of four order parameters when computing
We considered 138 collective variables consisting of oxygen–oxygen distances; oxygen–hydrogen distances for initially bound water molecules; all angles formed by Oλ and its four closest oxygen neighbors; and the Steinhardt order parameters of orders 3, 4, and 6 (32) (see Materials and Methods for more details). In addition, the order parameters already considered were added. Fig. 5A shows the resulting decision tree. Remarkably, of all of the input parameters, the
Results from the machine-learning analysis. (A) Classification and regression tree for predicting the outcome of initiated trajectories. Here, we considered several additional collective variables (description in Materials and Methods), but only a small subset is eventually needed for constructing the tree:
Conclusions
We investigated the autoionization of water at room temperature, using an unconstrained ab initio rare-event simulation method. Our simulations sample reactive events that happen on the timescale of minutes and we demonstrated that autoionization can be initiated by the hypercoordination of a stretched OH bond and the compression of a hydrogen bond wire as suggested by Hassalani et al. (17). However, these are not sufficient conditions. Only when the wire is strongly condensed (
Due to the multiple correlated factors that influence the water autoionization, we combined our analysis method with ML techniques which identified additional parameters not considered before, in particular the O–H stretch of the oxygen closest to Oλ. Even though the ML result did not outperform the level of predictiveness by the human effort based on intuition, visual inspection of many molecular movies, and intensive trial-and-error approaches, the ML approach found all previously identified parameters very efficiently and, in addition, revealed some equally important parameters that were overlooked. We therefore believe that ML applied to path sampling has a great potential especially since data limitations will become less of an issue in the future due to the further expected increase of high-performance computing, a better parallelization scheme of sampling unequal trajectory-length path ensembles, and the use of more efficient Monte Carlo (MC) path-generating moves (37). It would therefore be promising to apply the same method to other aqueous-phase chemistry studies which so far have mainly been based on biased dynamics (31, 38).
The fundamental understanding of reaction triggers that can be gathered by this approach could open up avenues of practical applications. For instance, even if not all identified parameters correlating with reactivity will necessarily imply causal correlation, it is plausible that an intelligent manipulation of their equilibrium distribution via external electric fields (39) or inclusion of additives might lead to catalytic ways to steer reactions and in particular water dissociation.
Materials and Methods
Simulation Methods.
The MD simulations required by the RETIS algorithm (14) were performed with the Born–Oppenheimer MD capabilities of the CP2K program package (40). We used the Becke–Lee–Yang–Parr (BLYP) functional with a DZVP-MOLOPT (41) basis set and a plane-wave cutoff of 280 Ry. The BLYP functional gives a reasonable description of the structure and dynamics of liquid water (42, 43) and the absence of dispersion corrections (44) is likely of minor importance for ion–water interactions where the dominant interactions are mainly electrostatic. However, we note that the BLYP functional is known to give an overstructured description of liquid water with a low diffusion coefficient (45). Previous studies on the recombination mechanism for water (17, 46) and for weak bases in water (18) have, however, found that the collective compression of the hydrogen bond wire and the motion of the protons are reproduced with different choices of the functional and basis set.
The initial system consisted of 32 water molecules placed in a cubic simulation box of
The transition region was divided into 20 path ensembles by positioning RETIS interfaces at
We performed 24,000 MC moves for each path ensemble, using the RETIS algorithm. This generated between 8,000 and 18,000 distinct trajectories in each path ensemble. The length of the trajectories ranged from 13.5 fs to 1,365 fs and we disregarded the first 400 trajectories in our analysis.
Analysis of Trajectories.
Crossing probabilities along the reaction coordinate λ were computed by matching the results of the different path ensembles. Projection of the crossing probability along
For trajectories harvested with the RETIS algorithm we calculated additional collective variables: the hydrogen bond wire length (
The orientation order parameter measures the distortion from a tetrahedral orientation of four water molecules around a central molecule and is defined by (27, 28)
After calculating these additional collective variables, we analyzed the trajectories using the methodology of van Erp et al. (19). For the analysis we used 100 subinterfaces for both
The classification models were constructed using CARTs (20) available within the R (49) software package. The mean of sensitivity and specificity was used as the classifier performance measure (50).
For the CART models we considered several sets of collective variables and we obtained these variables at the frame in the trajectories where the order parameter first crossed 1.15 Å. The trajectories were classified as reactive if they reached a
In the best-performing model (performance measure for training 0.89 and for testing 0.88) we considered 138 collective variables: all oxygen–hydrogen distances for initially bound water molecules, all oxygen–oxygen distances involving Oλ, the averaged distances between Oλ and its
Acknowledgments
The authors thank Øivind Wilhelmsen for fruitful discussions. The authors thank the Research Council of Norway Projects 237423 and 250875 and the Faculty of Natural Sciences and Technology, Norwegian University of Science and Technology (NTNU) for support. This research was supported in part with computational resources at NTNU provided by the Norwegian Metacenter for Computational Science (NOTUR), www.sigma2.no.
Footnotes
↵1M.M. and A.L. contributed equally to this work.
- ↵2To whom correspondence may be addressed. Email: titus.van.erp{at}ntnu.no or anders.lervik{at}ntnu.no.
↵3Deceased December 27, 2017.
Author contributions: A.L. and T.S.v.E. designed research; M.M., A.L., E.R., and T.S.v.E. performed research; M.M., A.L., V.V., and B.K.A. analyzed data; and M.M., A.L., E.R., V.V., and T.S.v.E. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission. P.L.G. is a guest editor invited by the Editorial Board.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1714070115/-/DCSupplemental.
Published under the PNAS license.
References
- ↵
- Eyring H,
- Henderson D
- Stillinger FH
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Trout BL,
- Parrinello M
- ↵
- Trout BL,
- Parrinello M
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Geissler PL,
- Dellago C,
- Chandler D,
- Hutter J,
- Parrinello M
- ↵
- Hassanali A,
- Prakash MK,
- Eshet H,
- Parrinello M
- ↵
- Cuny J,
- Hassanali AA
- ↵
- van Erp TS,
- Moqadam M,
- Riccardi E,
- Lervik A
- ↵
- Breiman L,
- Friedman J,
- Olshen R,
- Stone C
- ↵
- ↵
- Piccini G,
- Alessio M,
- Sauer J
- ↵
- ↵
- Dellago C,
- Bolhuis PG
- ↵
- van Erp TS,
- Bolhuis PG
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Steinhardt PJ,
- Nelson DR,
- Ronchetti M
- ↵
- ↵
- ↵
- ↵
- ↵
- Riccardi E,
- Dahlen O,
- van Erp TS
- ↵
- ↵
- ↵
- Hutter J,
- Iannuzzi M,
- Schiffmann F,
- VandeVondele J
- ↵
- ↵
- ↵
- ↵
- ↵
- Gillan MJ,
- Alfè D,
- Michaelides A
- ↵
- Hassanali A,
- Giberti F,
- Cuny J,
- Kühne TD,
- Parrinello M
- ↵
- ↵
- Cormen TH,
- Leiserson CE,
- Rivest RL,
- Stein C
- ↵
- R Core Team
- ↵
- Brodersen KH,
- Ong CS,
- Stephan KE,
- Buhmann JM
Citation Manager Formats
Article Classifications
- Physical Sciences
- Chemistry