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 H₂O species, λ is the largest covalent O–H bond distance, and when the system contains OH⁻ and H₃O⁺ species, λ is taken as the shortest distance between the oxygen in OH⁻ and the hydrogen atoms in H₃O⁺. In the following, we refer to the oxygen atom used for the order parameter as O^λ. The type of species (OH⁻, H²O, or H³O⁺) 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

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 H³O⁺ 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 |OH|1, |OH|2, |OH|3 are also indicated. These distances are used to investigate the possible concerted motion of hydrogen atoms along the wire. (C) The crossing probability (PA) and average energy of trajectories (⟨E⟩) as a function of the order parameter. The (black) dashed line is calculated using an alternative definition of the order parameter (λ′) where the trajectory length (in femtoseconds) defines the order parameter for λ>3 Å. The horizontal dotted-dashed line is the crossing probability (4.0×10−15) obtained for long paths (λ′≥1 ps). The activation energy is equal to the plateau value of the average energy which approaches 17.8 kcal/mol. The shaded area (1.15<λ<2.0) is the domain used for the predictive power analysis.

From our RETIS simulations, the water dissociation rate constant, kD, can be obtained as the product of a flux, fA, and a (conditional) probability, PA(λN|λ0):kD=fA×PA(λN|λ0).[1]Here, λ0 and λN are interfaces defining the initial (λ<λ0) and final (λ>λN) states and PA(λN|λ0) is the probability of reaching the final state before (possibly) reentering the initial state, given that the initial interface λ0 has been crossed. The flux, fA, is a measure of the frequency of crossings with λ0. Since we consider the dissociation of any water molecule in our system, we have normalized fA by the number of water molecules present. Typically, for rare events, the crossing probability is very small and in practice, PA(λN|λ0) is calculated by first positioning several more interfaces λ0<λ1<…<λN between the initial and the final state. The overall crossing probability is then obtained as a product of several (history-dependent) conditional probabilities (14). The conditional probabilities are calculated in a separate path ensemble simulation where the [i+] path ensemble defines the collection of paths crossing λi. The number and location of the interfaces alter the efficiency of the method, but not the results.

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 H²O 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, λ′, which equals the trajectory length (in femtoseconds).

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 H₃O⁺ 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 λ>3 Å, we consider λ′ as the order parameter and we have used a threshold of λ′≥1 ps as a criterion to identify a stable dissociation event. This choice is rather arbitrary since there is not a clear separation of timescales for the reverse recombination reaction which would result in another flat plateau region of the crossing probability. With a threshold of 1 ps, the crossing probability is PA=4.0×10−15. Combined with the initial flux, calculated to be fA=2.9×10−3 fs⁻¹ in our simulations, the resulting dissociation constant is kD=fA×PA=1.1×10−2 s⁻¹. An alternative rate constant not requiring any time threshold can be defined by counting the trajectories that undergo a hydrogen swap; i.e., in the last frame some of the water molecules have swapped their protons. The rationale behind this definition is that the proton swap must imply a significant reorganization of the hydrogen bond network so that the reverse reaction can be considered as an independent recombination reaction. Vice versa, the forward reaction has established a quasi-stable state since it is not followed up by a correlated reverse reaction. This definition yields a rate of kD=0.16 s⁻¹.

Comparing with experimentally determined dissociation constants at 25 °C [kD=2.5×10−5 s⁻¹ (7) and kD=2.04×10−5 s⁻¹ (8)] we overestimate the rate constant by a factor 500 (although the simulated rate will drop and gets closer to the experimental rate if a larger threshold is chosen). Considering all factors that play a role in the accuracy (statistical error, functional, small system size, purely classical treatment of protons, the time threshold value) the deviation with experiments is satisfactory and comparable to other density functional theory studies. Depending on the functionals considered in the ab initio calculation, energy barriers may be in error by 10–20 kJ/mol (22) which at room temperature would already correspond to a factor 55–3,000 difference between experimental and theoretical rate constants. Still, density functional theory generally manages to reproduce trends and mechanistic information in reasonable agreement with experiments (23).

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 ξ=(ξ1,ξ2,…). In principle, these ξis can be functions of all positions and momenta in the system, and they do not necessarily have simple physical interpretations. Since the ability to form hydrogen bonds is one of the characteristic features of water (26) and since previous computational studies have demonstrated the relevance of the hydrogen bond wire connecting the ionized species (16, 17), we have focused on a set of relatively simple collective variables which quantify the hydrogen bond network and the distortion from tetrahedral geometry.

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 i−1 other water species at the first point in time when λ is greater than a given threshold value, λt=1.15 Å. Typically, this threshold is reached within 3–6 fs in our trajectories. This defines a wire containing i water species whose length, wi, is obtained as the sum of the O–O distances of consecutive members.

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, na, and (ii) donated, nd, by the water species containing O^λ; (iii) the tetrahedral order parameter, q, obtained using the angles defined by O^λ and its four nearest oxygen atoms (27, 28) (by the definition q=1 for a perfect tetrahedral structure and q≠1 otherwise); and (iv) an angle order parameter, qcos, defined as the smallest of the cosine of the two internal angles in the wire. We refer to Materials and Methods for additional information on these collective variables.

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 λr and λc defined such that λr>λc>λ0. A trajectory is considered reactive if it reaches the specified λr; otherwise it is considered nonreactive. At the first crossing point with λc, we record the ξs and form two distributions using the reactive/nonreactive classification: rλc,λr(ξ), the fraction of λc-passing trajectories that cross λc at a point ξ and reach λr, and uλc,λr(ξ), the fraction of λc-passing trajectories that cross λc at a point ξ but fail to reach λr. These two distributions give information on the relation between the additional order parameters and the reactivity. For instance, if uλc,λr(ξ)=0, it could be that ξ is inaccessible, but if we can cross λc at ξ, the trajectory will be reactive. To quantify the importance of the different ξs, we calculate the predictive ability, TAλc,λr, defined as (19)TAλc,λr=1−1PA(λr|λc)∫rλc,λr(ξ)uλc,λr(ξ)rλc,λr(ξ)+uλc,λr(ξ)dξ,[2]such that 1≥TAλc,λr≥PA(λr|λc). If the collective variables do not correlate with reactivity, the lower limit is attained but if the ξs are relevant for the reaction, TAλc,λr>PA(λr|λc). We use the ratio TAλc,λr/PA(λr|λc)≥1 to measure how much the predictive ability is increased when considering the extra ξs, compared with using the crossing probability alone. Note that the definition in Eq. 2 shows that if the overlap of the two distributions is small, then the predictive ability increases.

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, w3, w4, w5) we find that w4 and w5 are more correlated with reactivity and that w4 is more relevant for larger λr (SI Appendix, Fig. S2). Thus, in the following we focus on wires containing four water molecules. For the water wires, we observe that when the ionic species are separated by at least two water molecules, the ionic state survives for a longer time compared with cases where they are separated by just one water molecule. This implies that (at least) three proton transfer events have occurred. We monitored the distances of the initially covalent O–H bonds and show these for the first (|OH|1), the second (|OH|2), and the third (|OH|3) transferred proton in Fig. 2. As can be expected from the Grotthuss mechanism (29, 30), the initial autoionization event is followed by several proton transfers in which the ionic species separate along the wire. Fig. 2 shows that this can happen in both a concerted and a stepwise way: The transfer of the first and the second proton occurs almost exclusively in a concerted way, while the transfer of the third proton (if it occurs) can happen stepwise or concertedly. This is also reflected in the waiting time between these events (SI Appendix, Fig. S3): The waiting time distribution between the second and the third proton transfer is broader compared with the first and the second transfer. To investigate the stability of the wires, we also calculated the hydrogen bond wire in time-reversed trajectories (SI Appendix, Fig. S4). We find that trajectories are indeed starting and ending with a contracted wire (w4<7.6 Å) as reported by Hassanali et al. (17), but at the end these wires do not necessarily contain the same oxygen atoms. This might occur due to an actual breakage of the hydrogen bond wire or by a lesser disruption (for example, by a shift of the selection of four consecutive oxygens within a five-membered wire). The majority of the longer trajectories reform via another wire, but there are still a significant number of long trajectories (>1 ps) for which the recombination is exactly the same as the dissociation path. This contradicts the hypothesis (16) that a breakage of the wire is a necessary condition to reach a metastable state. Also, visual inspection shows that relatively long trajectories exist in which the hydrogen bond wire remains intact except for some very short on/off fluctuations in the hydrogen bonds. We find the abovementioned hypothesis therefore difficult to defend. Conversely, we can also examine whether an actual breakage always leads to a long-lived metastable state. For this we adopt again the assumption that all trajectories with a hydrogen swap necessarily imply an indisputable breakage of the hydrogen bond wire. SI Appendix, Fig. S5, shows that trajectories with a proton swap are on average longer, but can still be relatively short (35 fs).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

The concerted behavior of the autoionization event. (A) The distances (|OH|i) of initially covalent O–H bonds for the first (i=1), the second (i=2), and the third (i=3) proton transfer in four trajectories. The arrows show the time direction and the different trajectories exemplify different types of hydrogen transfer: failed stepwise (dark-gray color, shown only for |OH|1–|OH|2), concerted (light-gray color), concerted only for |OH|1–|OH|2 (blue color), and concerted stepwise (orange color). (B) Using all trajectories in the final path ensemble, densities for |OH|1–|OH|2 and |OH|1–|OH|3 have been obtained (Left column). Right column shows the density when considering trajectories with a length tpath>60 fs. All trajectories were collected from the final path ensemble.

Comparing the additional collective variables (SI Appendix, Figs. S6 and S7), we find that nd is less relevant than the other variables and we do not consider it further. The other collective variables are more correlated with reactivity and in Fig. 3A, we show the predictive ability for some of their combinations. In Fig. 3B we show TAλc,λr as function of λr≤2 Å for λc=1.16 Å compared with the crossing probability using several combinations of the collective variables. Fig. 3 A and B shows that we can increase the predictive ability by a factor 107 compared with the crossing probability. We note that since the crossing probability is small in this case, with a TAλc,λr∼0.4, we cannot perfectly predict the outcome. This indicates that there are other collective variables important for the description, possibly even nonlocal ones, as suggested by Geissler et al. (16). Also, we stress that here we are focusing only on the first concerted-jump step of the reaction in which the order parameter increases from 1.16 Å to 2.0 Å. As is clear from Fig. 1C the vast majority of trajectories reaching λ=2.0 Å will not lead to long-lived metastable states. Predicting this from the very first snapshot seems still a step too far since it depends on collisions between water molecules after many MD steps far away from the initially stretched OH bond.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Increasing the predictive power for water autoionization by considering additional collective variables. (A) The predictive power (TAλc,λr[ξ]) relative to the crossing probability (PA(λr|λc)) using additional collective variables: hydrogen bond wire length (ξ=w4), the orientation order parameter (ξ=q), the angular order parameter (ξ=qcos), and the number of hydrogen bonds accepted (ξ=na) by the O^λ species. (B) The predictive power and the crossing probability as a function of λr for λc=1.16 Å and different combinations of collective variables. Due to the threshold criterion for defining the wires (main text), the probability is shifted so that PA=1 for λ<1.15 Å.

Inspecting the initiation conditions in more detail, we investigate the reactive and nonreactive distributions rλc,λr(ξ) and uλc,λr(ξ) in Fig. 4 for λc=1.16 Å and λr=2.0 Å. Here, we examine all dissociation events, even the ones that recombine quickly and show the distributions for ξ=(w4,na) in Fig. 4A [see SI Appendix, Fig. S8 for the distributions for ξ=(w4,q) and ξ=(w4,qcos)]. Along the w4 coordinate we observe a clear separation of the two distributions which indicates that trajectories crossing λc=1.16 Å have a larger probability of being reactive for shorter wires (smaller w4). This supports the hypothesis of a “compressed” wire as an important condition for autoionization, as first suggested by Hassanali et al. (17). Along the na coordinate we observe a higher probability of reactivity for wires in which O^λ is hypercoordinated, which was also proposed by Hassanali et al. (17). Still, the chance of not being reactive is larger at any point (w4,na) in Fig. 4A [rλc,λr(ξ) would not be visible if it had not been normalized]. For example, if (i) 7.15<w4<7.6 and at the same time na=3, the probability for a reactive event is 3.6⋅10−6, which is small but still a factor 58 larger than the chance to be reactive from a random point at λc. In a more extreme case, if (ii) w4<7.3 and simultaneously na=4, the chance increases to 0.15. The predictive ability TAλc,λr provides a weighted average of these chances in which the weights are proportional to the relevance (19); since of all reactive trajectories, 45% cross λc in region i and only 0.6% in region ii, the latter will have 75 times lower weight.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Initiation conditions and local collective variables. (A) Reactive (rλc,λr(ξ)) and nonreactive (uλc,λr(ξ)) distributions for ξ={w4,na} and λc=1.16 Å and λr=2.0 Å. For visualization purposes, the depicted distributions are normalized [implying magnification of 107 for rλc,λr(ξ)]. Top and Right Insets show the one-dimensional projections of the distributions. A clear separation of the two distributions can be seen along the w4 coordinate, indicating that reactive trajectories are more compressed compared with nonreactive trajectories. In addition, the oxygen atom used in the order parameter calculation (O^λ) accepts on average a larger number of hydrogen bonds in reactive trajectories, compared with nonreactive trajectories. (B) Illustrative snapshot from a reactive trajectory where O^λ is shown in blue. The four surrounding oxygen atoms which are used for the calculation of the tetrahedral order parameter q are shown in orange. The water wire is highlighted with a yellow line (and gray transparent spheres) and the angle parameter qcos is indicated. In this snapshot, the water wire is compressed, q exhibits deviation from a tetrahedral structure, qcos indicates that three oxygen atoms are lining up in the wire, and O^λ accepts three hydrogen bonds and donates one (shown with green lines).

If we consider the q coordinate, we observe that rλc,λr is shifted toward lower q values compared with uλc,λr, which indicates that a distortion from a tetrahedral arrangement around the dissociating water species may also initiate the event. This finding is somewhat surprising as in some other aqueous phase chemical reactions the opposite effect was found (31). Similar conclusions can be drawn for the distribution of ξ=(w4,qcos). Here, there is a peak along the qcos coordinate for the reactive distribution closer to a linear arrangement of the water molecules. In Fig. 4B we show a representative snapshot, obtained early (after 3 fs) in a reactive trajectory. Overall the results shown in Fig. 3 report that a compression of the water wire (measured by w4) and hypercoordination (measured by na) or distortion (measured by q and qcos) are necessary initiation conditions for autoionization. However, these are not sufficient conditions as shown by the values of TAλc,λr in Fig. 3B: Still 60% of the trajectories starting off within the ideal ξ parameter range fail to establish a concerted proton jump.

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 λ>1.15 Å is first attained, to predict the outcome. We used several ML techniques in which every odd path ensemble was included in the calibration and the even path ensembles were used for the test set. An alternative split in which the data within each path ensemble were evenly divided in two gave similar results. Moreover, as heavily skewed distributions are difficult to treat with ML, we further omitted the reweighting of the datasets with the statistical weights of the corresponding path ensembles. However, we applied the ML techniques as a qualitative approach to find new parameters that could be tested quantitatively within the predictive power method (19).

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 TAλc,λr. For instance, excellent predictive performances (>90%) were obtained using the ensemble-based gradient-boosting machines (35, 36). However, the interpretation of the model is problematic since an ensemble of 100–150 deep decision trees (added in a sequence) is used. Although the performance is improved, the chance of overfitting with accidental correlations increases. We have therefore restricted ourselves to the single-tree–based decision models based on classification and regression decision trees (CART) (20). The restriction to four order parameters for the TAλc,λr function is based on similar reasons. Adding more parameters gives more sparse matrices representing the reactive/nonreactive distributions, and, as a result, numerical integration for computing the overlap between these distributions becomes very sensitive to the bin size and could underestimate the overlap due to bins being empty by insufficient statistics.

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 w4 parameter is both on top of the decision tree and the most important variable as measured by the reduction in the classification error attributed to each variable at each split in the decision tree (20) (SI Appendix, Fig. S9). Also the tetrahedral ordering and the number of accepted hydrogen bonds appear in the decision tree. To describe the first effect, the ML approach prioritized the Steinhardt q4 order parameter above the similar q parameter previously used by us. Some distances that also appear in the decision tree like d25, the distance between O^λ and its 25th closest oxygen, are most likely due to accidental correlations caused by the limited size of the dataset. This is verified by inspecting the importance of this variable: d25 does not appear among the 20 most important variables (SI Appendix, Fig. S9), and, in fact, other similar variables (e.g., d24) are ranked higher, albeit with low importance. A more important and intuitively sound parameter that is suggested by the ML approach is λ2, the OH distance between the oxygen closest to O^λ and its hydrogen with the largest intramolecular bond. Recomputing the predictive ability using parameters from the ML tree (Fig. 5B) did not yield higher performances than the combination w4, q, na, and qcos, but should be conceived as equally good, considering statistical uncertainties.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

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: w4, q4 [the Steinhardt fourth-order parameter (32)], λ2 (the length of the stretched hydrogen bond in the water molecule closest to the O^λ species), di (the distance from O^λ to the ith closest oxygen), and d¯i (the average distance considering the i closest oxygens). The notation for the nodes is explained with the stand-alone node in the top left corner. This tree predicts trajectories to be reactive, i.e., reaching a λ≥2, or nonreactive based on the collective variables obtained at the frame in the trajectories when λ is first ≥1.15. The nodes predicting reactive trajectories are colored blue (class 1) while the nodes predicting nonreactive trajectories are colored green (class 0). Note that the percentages at the bottom of the squares do not reflect the physically correct fractions since path ensembles were not reweighted using their statistical weights. The rules are textual representations of traversing the tree; for instance, rule 5 (which predicts reactive trajectories) can be expressed as w4≥7.6 and λ2≥1.1. These rules give different initiation conditions, and they are listed in SI Appendix, Table S1, for the bottom row of nodes. (B) The predictive power and the crossing probability as a function of λr for λc=1.16 Å and different combinations of collective variables. Here we compare the predictive power using collective variables we identified with variables marked as important by the machine-learning analysis. (C) Reactive (rλc,λr(ξ)) and nonreactive (uλc,λr(ξ)) distributions for ξ={λ2,d¯2} and λc=1.16 Å and λr=2.0 Å. For visualization purposes, the depicted distributions are normalized. Top and Right Insets show the one-dimensional projections of the distributions.