Autoionization in Liquid Water

See allHide authors and affiliations

Science  16 Mar 2001:
Vol. 291, Issue 5511, pp. 2121-2124
DOI: 10.1126/science.1056991


The dissociation of a water molecule in liquid water is the fundamental event in acid-base chemistry, determining the pH of water. Because of the short time scales and microscopic length scales involved, the dynamics of this autoionization have not been directly probed by experiment. Here, the autoionization mechanism is revealed by sampling and analyzing ab initio molecular dynamics trajectories. We identify the rare fluctuations in solvation energies that destabilize an oxygen-hydrogen bond. Through the transfer of protons along a hydrogen bond “wire,” the nascent ions separate by three or more neighbors. If the hydrogen bond wire connecting the two ions is subsequently broken, a metastable charge-separated state is visited. The ions may then diffuse to large separations. If, however, the hydrogen bond wire remains unbroken, the ions recombine rapidly. Because of their concomitant large electric fields, the transient ionic species produced in this case may provide an experimentally detectable signal of the dynamics we report.

A randomly chosen, intact water molecule will dissociate in liquid water, producing hydronium (H3O+) and hydroxide (OH) ions, within ∼10 hours (1). The event is thus extremely rare on the femtosecond time scale of molecular motions. But when it occurs, the system crosses a transition state region quickly, generating ions that are separated by ∼6 Å (2). Long ago, Eigen and De Maeyer argued that from this intermediate state, diffusion of the ions away from each other is relatively facile (1). The dynamical bottleneck to autoionization is therefore associated with the initial formation of this intermediate state. Eigen and De Maeyer further conjectured a microscopic mechanism by which this state is accessed. In doing so, they recognized that the important coordinates for chemistry in liquid water differ from those appropriate for gas phase reactions. Specifically, they imagined that dissociation of a water molecule is driven by rearrangement of solvating water molecules. In this picture, the transition state involves specific solvation structures and hydrogen bond patterns that accommodate hydronium and hydroxide ions. But even today, the microscopic motions of the system as it crosses such a transition state, on the subpicosecond time scale consistent with Eigen and De Maeyer's kinetic model, cannot be resolved experimentally. Until now, therefore, the role of solvent fluctuations in autoionization has not been demonstrated, nor have its molecular details been determined.

Molecular dynamics simulations would seem well-suited to explore atomic motions on the time scale of transition state crossing. With current computers, however, reasonably accurate modeling of chemical dynamics in condensed phases can generate trajectories lasting only picoseconds, far short of the typical times between dissociation events. Straightforward dynamics simulations thus cannot bridge the huge separation of time scales. An attempt to avoid this problem by controlling a proposed reaction coordinate, specifically the distance between a particular oxygen-hydrogen pair, failed to produce stable separated ions (3, 4). This artificially chosen coordinate, which is appropriate for dissociation in the gas phase, is evidently not the true reaction coordinate for dissociation in liquid water.

We show here that the fluctuations driving autoionization involve the concerted motion of many atoms, as is expected for a condensed phase chemical reaction. This insight is obtained by using Car and Parrinello's molecular dynamics (CPMD) (5) and transition path sampling (6–9). Although charges may be separated artificially in water by moving a small number of atoms, the driving force for autoionization can only be described with collective coordinates. These coordinates, as determined from our simulations, are not of the specific nature that Eigen and De Maeyer imagined. In particular, they do not describe unique arrangements of solvent molecules. They are instead more akin to coordinates that appear in Marcus's theory of electron transfer (10), describing, for example, the electric field generated by solvent molecules. An ensemble of autoionization transition states thus exists, characterized by certain values of collective solvent coordinates but diverse in the positions of individual atoms.

Combining CPMD with transition path sampling enables the separation of time scales in autoionization to be overcome. Based on electron density functional theory, CPMD allows efficient classical nuclear dynamics simulations in condensed phases with forces determined at each time step by quantum chemistry calculations. From comparisons of force fields relevant for proton transfer (11) and dissociation (12) in water clusters, one expects CPMD [using the BLYP functional (13, 14)] to be reasonably accurate for autoionization in liquid water.

The nuclear dynamics that we report are classical, but the coordinates we identify as important for autoionization have low associated masses, and quantum mechanics may thus influence kinetics strongly. In particular, electric field fluctuations in water arise in large part from high-frequency oscillation of local dipoles (i.e., librational motions). Quantization of these motions can substantially enhance rate constants, as in the case of electron transfer (15). Nevertheless, just as in electron transfer, these effects of tunneling and zero-point motion should not change the qualitative nature of the reaction coordinate. Indeed, for transport of an excess proton in water, it has been observed that no fundamentally different pathways are made accessible by tunneling of a few light degrees of freedom (16). We therefore expect that our classical nuclear trajectories explore the regions of configuration space pertinent for autoionization.

In the trajectories that we obtained, the simulated system consists of 32 water molecules periodically replicated in space. The dimensions of our central simulation cell are therefore not large, and we cannot probe the influence of correlations beyond a length of ∼10 Å. Long-range forces were not, however, truncated artificially. Instead, interactions among all periodic images were computed with Ewald summation. Calculations of the structural properties of liquid water at equilibrium using a nearly identical model system (17) suggest that the most important long-range fluctuations are captured in this representation. More importantly, a system of this size is sufficient to describe the dissociation thermodynamics of a pair of simple, monovalent ions in water (18).

Transition path sampling, with which our trajectories are gathered, is an importance sampling in trajectory space. Attention is restricted to trajectories that cross the transition state surface on a molecular time scale (here, ∼200 fs). This is accomplished through a Monte Carlo sampling of the ensemble of transition paths. From a given reactive pathway, we generate a trial pathway by a small displacement, such as changing the atomic momenta slightly at some point along an existing path. A new path is then obtained by integrating the equations of motion forward and backward in time. The resulting trajectory is accepted if it still leads from reactants (here, neutral water molecules) to products (here, separated ions); otherwise, it is rejected. Only small momentum displacements are required because the chaotic nature of dynamics ensures that the resulting trial trajectory will differ substantially from the original path. Indeed, only a few such moves are typically needed to produce a statistically independent trajectory (19). Our sampling of autoionization pathways thus represents ∼10 uncorrelated reactive events, harvested with their physically appropriate weights.

In order to sample transition paths, it is necessary only to distinguish between reactants and products and to generate an initial reactive pathway. Unlike conventional studies of barrier crossing dynamics, transition path sampling does not require the specification of a reaction coordinate. This feature of the method is crucial for studying complex systems such as water, in which no simple coordinate describes the motion of the system through an ensemble of different transition states.

In a typical trajectory from our sampling of transition paths (Fig. 1), the system travels from the neutral state, through a transition state, into the charge-separated state within ∼150 fs. Basins of attraction and transition states must be identified statistically in this high-dimensional system. In particular, from any configuration, we may initiate many trajectories, with initial momenta chosen at random from an equilibrium distribution. If nearly all of these trajectories lead to neutral water molecules, then the configuration belongs to the neutral state. Similarly, if most trajectories lead to separated ions, then the configuration belongs to the charge-separated state. A configuration is considered to be a transition state if equal numbers of trajectories result in neutrality and charge separation (20, 21). On the basis of this criterion, the transition state surface for autoionization lies between the configurations shown in Fig. 1, D and E.

Figure 1

(A through F) Snapshots from an autoionization trajectory for 32 water molecules periodically replicated in space. The illustrated configurations are separated by 30 fs. Dynamics were computed with kinetic energy consistent with room temperature and with a force field determined by the BLYP electronic density functional theory. The trajectory was generated by the CPMD algorithm and was harvested by transition path sampling. Red spheres represent oxygen atoms, and white spheres represent hydrogen atoms. The yellow and blue oxygen atoms highlight the hydronium and hydroxide ions, respectively. The dotted lines show hydrogen bond wires along which ionic species move relatively quickly.

The exhibited trajectory begins in the neutral state (Fig. 1A). As in a typical equilibrium configuration of water, each oxygen atom is chemically bound to two hydrogen atoms. Dissociation begins when a fluctuation in solvent electric field causes the cleavage of an oxygen-hydrogen bond on the nascent hydroxide ion. The unbound proton is transferred to a neighboring water molecule that becomes, transiently, the hydronium ion. Within 30 fs, the hydronium ion shifts by proton transfer, so that the two ions are no longer nearest neighbors (Fig. 1B). This transport of ions by proton transfer rather than molecular translation, known as the Grotthuss mechanism (22,23), is essential for autoionization. Over the next 60 fs (Fig. 1, C and D), the solvent-generated electric field continues to push the ions away from each other. Both ions move rapidly along a path of hydrogen bonds, effectively a wire for the conduction of protons. In the following 30 fs, a crucial fluctuation carries the system over the transition state; namely, the hydrogen bond wire that connected the ions is broken through reorganization of the hydrogen bond network (Fig. 1E). Disconnected, the ions can no longer recombine by rapid proton transfer along existing hydrogen bonds. After this event, the system resides in the charge-separated state (Fig. 1F). The ions are accommodated by the aqueous environment, as their characteristic solvation structures (24–26) are formed.

Distinguishing between the neutral reactant state and the charge-separated product state for autoionization involves more than just the interionic distance. In our simulated system, hydronium and hydroxide ions can be separated by a maximum of ∼10 Å. At such microscopic separations, it is possible that the system remains in the basin of attraction of the neutral state. For example, recombination occurs within 100 fs after we artificially remove a proton from one water molecule in a neutral system and place it on a second molecule several molecular diameters away. The fast recombination occurs by proton transfer along a hydrogen bond wire. We observed this phenomenon for several pairs of water molecules. The interionic distance is therefore not a useful order parameter for distinguishing the neutral and charge-separated states for autoionization.

When a strong potential difference exists along a hydrogen bond wire, as in the case of separated ions, proton motion can occur without a substantial barrier. When no short hydrogen bond wire connects the ions, however, conduction of protons between them cannot occur by small-amplitude motions. Recombination then requires rearrangement of the hydrogen bond network. Indeed, rapid recombination does not occur when a proton from one water molecule is artificially placed on another while simultaneously ensuring that the resulting ions are not connected by a short hydrogen bond wire. Several trajectories initiated in this way remained in the charge-separated state for longer than 1 ps before a hydrogen bond wire formed and the ions recombined. (One of these trajectories, reversed in time, was used to begin our path sampling.) This picosecond lifetime is consistent with the rate at which a hydrogen bond in liquid water changes orientation from one water molecule to another (27–29). In the metastable charge-separated state, the ionic structures move by proton transfer and reorganization of hydrogen bonds. But the ions do not move systematically toward each other. In an extended system, they may in fact diffuse to large separations, on the picosecond time scale of the Grotthuss mechanism (24–26, 30, 31). In our periodic system, however, a short hydrogen bond wire eventually forms, and rapid proton transfer neutralizes the system.

To describe the fluctuations that control available routes for proton transfer, we define an order parameter l as the number of bonds in the shortest hydrogen bond wire connecting the ions. In our finite simulation cell, l is effectively constrained to be less than 5. Longer hydrogen bond wires are undoubtedly present in bulk water and may even contribute to the transition state ensemble, influencing the kinetics of autoionization. In our simulated system,l is most often between 3 and 5 in the charge-separated state, changing as the ionic structures move and as hydrogen bonds are broken and formed. If l = 2, recombination usually occurs within 100 fs. The reactant and product states appear to be characterized by l = 0 and l ≥ 3, respectively. These criteria, which we use to distinguish between the two states in transition path sampling, are consistent with the critical ion separation of 6 Å estimated by experiment (1). This length coincides with a typical spatial distance spanned by a hydrogen bond wire of length l = 3. In fact, there is a distribution of hydrogen bond wire lengths at which the ions can become stably separated. For each length, diffusion of the ions away from each other will occur with a different rate constant. A careful examination of weak acid dissociation kinetics should therefore reveal complex nonexponential behavior.

Although the length of the hydrogen bond wire distinguishes between neutral and charge-separated states, it does not fully describe the reaction coordinate. Some configurations with l = 2 belong to the reactant basin of attraction, whereas others belong to products. A second coordinate, related to solvation fluctuations, is needed to distinguish between these groups of configurations and to specify a transition state. The nature of this coordinate is exhibited by computing the potential of protons in the hydrogen bond wire, with all other degrees of freedom held fixed at a particular configuration. When an appropriate solvent electric field is present, this potential has two local minima, corresponding to separated ions and neutral molecules. We located these minima using molecular dynamics (32). We then probed the energy of configuration space between them by interpolating linearly between their structures (Fig. 2). The energy difference between the two minima, ΔE, indicates the solvent preference for separated ions over neutral molecules, effectively the direction and strength of the solvent-generated electric field.

Figure 2

(A through C) Potential of protons in the hydrogen bond wire connecting hydronium and hydroxide ions. A one-dimensional projection of the 3 ×l–dimensional potential is obtained by interpolating between the structures of the two local minima. These minima were determined by quenching wire protons at low temperature while holding all other coordinates fixed. Access to different wells of the potential was provided by moving wire protons artificially, so that ions or neutral molecules were formed. The coordinate q is defined to interpolate between the ionic and neutral configurations,r(q) = (1 –q)r (ion) +q r (neut), wherer (ion) and r (neut) are the 32 × 3 × 3–dimensional position vectors specifying the configurations of the entire system in the ionic and neutral minima, respectively. We define the coordinate ΔE as the difference in energy between the two minima, i.e., ΔE= E[r(1)] – E[r(0)], where E[r(q)] is the energy of configuration r(q). In (B) and (C), the left minimum (at q = 0) corresponds to separated ions. The right minimum (at q = 1) corresponds to neutral water molecules. These ionic and neutral structures are labeled in (A), (B), and (C) by schematic wire configurations, with colors as in Fig. 1. The potentials in (B) and (C) were calculated for configurations withl = 2, for which the solvent electric field favors separated ions. This electric field is not present in the configuration corresponding to (A). Because the ions are unstable in this case, r (ion) cannot be obtained by quenching. Instead, it was constructed artificially by moving wire protons to give reasonable ionic structures. All energies are plotted relative to the energy of the neutral minimum,E neut.

The proton potential so obtained is plotted in Fig. 2 for three configurations along an autoionization trajectory. The first configuration (Fig. 2A) belongs to the neutral state. Because the bond-destabilizing electric field has not yet appeared, separated ions are not stable. A single minimum is present, corresponding to neutral water molecules. When the electric field begins to appear (Fig. 2B), the separated ions become metastable with respect to proton motion. Here, the ions are just 2 kcal/mol higher in energy than neutral molecules. As the electric field fluctuation increases, the separated ions become still more energetically stable, at the expense of neutral molecules. In Fig. 2C, the ions are more stable by over 20 kcal/mol. The solvent has not only lowered the barrier to bond dissociation, but also has driven protons along a hydrogen bond wire to form hydronium and hydroxide ions.

We describe the solvent coordinate that stabilizes ions, ΔE, as an electric field because it arises primarily from long-range electrostatic interactions. Local properties that we examined, such as ion coordination number and the presence of specific hydrogen bonds, fail to account for the bond-destabilizing fluctuation in our simulations. Furthermore, the stabilization of ions indicated inFig. 2, B and C, is diminished substantially when we artificially remove outer coordination shells of the ions. Thus, ΔEdoes not arise from a few nearby water molecules. Instead, it is analogous to the collective coordinates that others have imagined for electron transfer (10) and for acid-base proton transfer (33–36). As in Marcus's theory of electron transfer, it is a rare solvent fluctuation that drives the motion of charges. In detail, however, ΔE differs from the previously defined coordinates, namely the solvent polarization field and the energy gap between diabatic bonding states. ΔEinvolves only the energy required to transfer protons adiabatically. The second component of the reaction coordinate that we identified, the hydrogen bond wire length l, is also analogous to a coordinate in these theories, namely the distance between ions. But in the case of autoionization, the appropriate separation coordinate describes the hydrogen bond wires that link the ions, rather than simply describing the distance between them, emphasizing the importance of connectivity in the hydrogen bond network.

We conclude that the dynamics of both electric fields and hydrogen bonding play important roles in the autoionization mechanism. Rare electric field fluctuations drive the dissociation of oxygen-hydrogen bonds. Ions produced in this way usually recombine quickly because the solvation fluctuation vanishes within tens of femtoseconds. But when such a fluctuation is coincident with breaking of the hydrogen bond wire (a process normally occurring about once every picosecond), rapid recombination is then not possible. It is with this coincidence of events that the system crosses a transition state. This scenario implies the existence of many short-lived hydronium and hydroxide ions in water. Decay of this transient population over ∼100 fs is an interesting and, in principle, observable signature of the dynamics revealed by our simulations.

  • * To whom correspondence should be addressed. E-mail: chandler{at}

  • Present address: Institute of Organic Chemistry, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland.


View Abstract

Navigate This Article