An analysis of hydrated proton diffusion in ab initio molecular dynamics

The pursuit to understand the mechanistic details of excess proton transport in water dates back to at least over 200 years ago. In an 1806 paper,¹ von Grotthuss first proposed a mechanism in which proton transfer can occur by breaking and forming chemical bonds of proximal water molecules so that the protonic defect diffuses or hops through the hydrogen-bond (H-bond) network. His contribution is remarkable considering that atomic theory was not well-understood at the time and as widely accepted as it is today. In addition to proton hopping between water molecules, the hydrated proton may also be transported with a water molecule as a hydronium cation (H₃O⁺) for extended periods of time giving rise to the more commonly known vehicular transport mechanism. Vehicular transport and the proton hopping mechanism, the latter now known as the Grotthuss mechanism,² play pivotal roles in transporting protonic charges in a vast array of aqueous, material, and biological systems. As such, it is critically important to understand the fundamental processes governing the inherently multiscale mechanism whereby dynamic chemical bond formation is coupled to the surrounding H-bond network, which are both coupled to larger-scale, slower motions in the system of interest.

As a community, our understanding of mechanistic details governing proton (and similarly hydroxide) transport has advanced tremendously in the last 25 years with the development of high-resolution spectroscopic measurements³ and computer simulations (some of the many examples^4–9 and reviews^10–12 that summarize our understanding of proton transport by simulation are available in the references). A reasonable model for studying proton transfer between water molecules must be able to account for the constantly changing covalent bonding topology in the H-bond network. The molecular simulation of such complex phenomena requires computational methodologies that are accurate to obtain robust estimates of physical observables and establish connections with experimental characterizations. For the case of proton transport in aqueous systems, the ab initio molecular dynamics (AIMD) simulation method is an instinctive choice for modeling the complex, dynamic rearrangement of covalent bonds and solvation structures. However, the computational cost for the explicit accounting of electronic degrees of freedom and inherent approximations in choice of electronic structure method (e.g., the exchange correlation functional and inclusion of van der Waals interactions in density functional theory (DFT)) can plague simulation results by requiring small systems (finite size effects), short simulation times (unconverged properties and/or poor equilibration), and known (and unknown) artifacts introduced by the level of theory employed. It is thus critical to sample sufficiently long AIMD trajectories to understand the robustness of statistical quantities related to events over the time scales of typical AIMD simulations (tens of picoseconds). To access larger scale phenomena, computationally efficient models that incorporate the essential physics and permit the dynamic evolution of the chemical bonding topology must currently be employed at varying levels of approximations (e.g., in no particular order, multistate empirical valence bond (MS-EVB),^7,13–16 multiconfigurational,¹⁷ ReaxFF,^18,19 Q-hop,²⁰ or SCC-DFTB).^21–23 

Despite the impression that an explicit inclusion of electronic degrees of freedom as in AIMD may yield simulation results in excellent agreement with experimental observables, that is generally not the case for bulk water systems. Commonly employed DFT water models, defined by the employed exchange-correlation potential, are found to have significant issues at room temperature. For the case of the Becke-Lee-Yang-Parr (BLYP) functional (with and without dispersion corrections), the computed melting temperature was found to be significantly higher than 300 K.^24,25 Consistent with a supercooled liquid, classical AIMD simulations with common DFT models result in an H₂O liquid that is over-structured and self-diffusion constants smaller than the experimental value (0.23 Ų/ps) by up to an order of magnitude at 300 K (∼0.01 Ų/ps), even though the Hamprecht-Cohen-Tozer-Handy (HCTH/120) functional yields a diffusion constant nearly half that of the experimental value.²⁶ Reported values for the self-diffusion coefficient for several other DFT approximations can be found elsewhere.^26–29 AIMD simulations of water at a higher temperature (315–330 K) have been found to lead to improved agreement in radial distribution functions (RDFs or g(r)) and diffusivity compared to experimental measurements at 300 K, but simulating a different thermodynamic state relative to the experimental conditions of interest is difficult to justify and can lead to unexpected side-effects and misinterpretation of findings, particularly as the system being investigated increases in complexity. This is not to imply that AIMD simulations have not contributed invaluable knowledge to the structure and dynamics of aqueous systems, but it is important to recognize these faults and be aware in a general sense of the pros and cons of simulating a particular system with a specific level of theory (e.g., DFT). For example, even though the known slow water dynamics will renormalize the time scales of the vehicular and hopping mechanisms of proton diffusion, because AIMD simulations can naturally account for charge transfer and electronic polarization, they can still provide useful insights into the mechanism of proton transport and its coupling to the surrounding H-bond network (the Grotthuss mechanism). Such detailed and statistically meaningful comparisons become even more critical with the recent implementation of periodic post-Hartree-Fock methods³⁰ that are now being used to examine bulk liquid systems.

There is similarly a large body of work on modeling proton transport in aqueous systems using AIMD and reactive simulation methods.^{4,5,10,12,26,28,31–33} In previous AIMD simulations, the importance of sequential H-bond cleavage for water molecules accepting and donating the excess proton, respectively, was observed as conditions before successful transfer of the proton with subsequent formation of a H-bond with the newly formed water molecule.²⁸ The present study is one example of how improved statistical sampling can provide new insight and refine a mechanistic picture for such a complicated process, like proton transport. In work published several years ago, it was shown from MS-EVB simulations that the proton transport in water is a process depending on the collective rearrangement of several water solvation shells.⁸ This picture was subsequently confirmed experimentally.³⁴ In the recent literature, Hassanali et al.^31,32 using AIMD simulations have observed the formation of and rapid proton transport along compressed water wires, periods of burst-and-rest dynamics, and directional H-bond and umbrella mode correlations with proton transfer events. The directionality of the surrounding H-bond network likely plays a critical role in the transport of protonic defects and water molecules.^8,31 H-bond rearrangements observed in these simulations have similarly important consequences in aqueous clusters,³⁵ ice,^36,37 interfaces,^38,39 and confined environments.^40–42 

It is well-documented that an excess hydrated proton can be trapped in its local solvation structure on the picosecond time scale as determined from both AIMD and MS-EVB simulation studies.^9,15,31,32 For reasons that were not previously clearly stated, the excess proton can “burst” and hop a distance of several water molecules on a much shorter time scale before entering another “rest” state. In order to build more physically accurate and refined computationally efficient models, it is important to understand and identify the key molecular features that such models should embody, whether they be derived purely from molecular mechanics Hamiltonians or constructed as empirically corrected density functional approximations. Building upon the large body of AIMD work to develop such models, it is also important to validate the statistical weight of transient processes correlated with both successful and unsuccessful proton transfer events and to also understand what effect varying simulation parameters may have on such computed physical observables.

The main goal of the present work is to understand the mechanisms governing the observed burst-and-rest behavior and their strength of coupling to the surrounding H-bond network as a function of both density functional approximation (BLYP, BLYP–D3, HCTH, B3LYP) and temperature (300 and 330 K). By repeating all simulations at two different temperatures, one can show how the increase in temperature changes the structural and dynamical properties in different AIMD setups and more easily establish connections with previous literature, which commonly simulates aqueous systems at elevated temperatures in an attempt to avoid the unphysical features of AIMD simulations which are based on most common DFT approximations. With connections to previous AIMD simulations, analysis of the hydrogen-bond connectivity is used to identify the role of two key features observed in proton transport simulations: concerted hops and weak hydrogen-bond donors. The question of concerted hops seeks to identify those time scales over which long-range proton transport can be classified as a series of sequential water hopping events and whether they always lead to long-range proton migration. The coupling of the fourth neighboring water molecules solvating the protonated water molecule (hydronium) before and after proton transfer events, which are weak H-bond donors, with the observed burst-and-rest proton dynamics is a key property investigated. Periods of burst-and-rest activities related to a series of hops onto nearby water molecules that may or may not be already pre-solvated by this fourth water (consistent with pre-solvation concepts)^10,28,43 are analyzed to identify the relative strength of correlations with the surrounding H-bond network, which ideally are physical observables that accurate and computationally efficient models should recover.

The Born-Oppenheimer AIMD simulations examined in this work were calculated using the Quickstep module in the CP2K software package.⁴⁴ As indicated in Table I, four density functional approximations were considered: BLYP,^45,46 BLYP with Grimme dispersion corrections (herein denoted BLYP–D3),^47,48 HCTH/120,⁴⁹ and the hybrid functional B3LYP.^50,51 The Gaussian and plane waves (GPW) method⁵² was employed using the Gaussian DZVP and TZV2P basis sets and plane wave cutoffs of 300 and 400 Ry for the auxiliary density. Goedecker-Teter-Hutter pseudopotentials⁵³ were used for the core electron states. All classical AIMD simulations contained one excess proton in 128 water (H₂O) molecules in a simple cubic simulation box of length 15.66 Å (1.0 g/cm³ density). The orbital transformation (OT) method was used to optimize the wave function at each step of the simulations with a convergence criteria of 1 × 10⁻⁶ a.u. The time step for integration was 0.5 fs. Configurations from previous AIMD simulations¹⁶ were used to initiate new simulations (Table I) in the constant NVT (canonical) ensemble using the canonical sampling through velocity rescaling (CSVR) thermostat with a relaxation time of 1 ps.⁵⁴ The DZVP–300Ry settings with the CSVR thermostat were selected to more readily facilitate comparison with recent AIMD studies.^31,32 For the BLYP–D3–TZV2P simulations, a Nosé-Hoover chain thermostat with characteristic frequency of 3000 cm^-1 was coupled to all degrees of freedom as used in previous studies.¹⁶ For all simulation settings with non-hybrid functionals, two target temperatures were selected for the NVT simulations, namely, 300 and 330 K. The final state from all NVT simulations were subsequently used to initialize simulations in the constant NVE (microcanonical) ensemble. The final states of the NVE simulations with BLYP–D3–DZVP were used as the initial states for the NVT simulations with the hybrid B3LYP–D3 functional, which was equilibrated at a single target temperature of 300 K. For simulations with the hybrid B3LYP–D3 functional, the libint library⁵⁵ was utilized and the settings of DZVP–300Ry for the basis set and truncation of Hartree-Fock exchange at 6.0 Å were selected as a balance between convergence and computational efficiency.^56,57 The accumulated simulation time of all AIMD trajectories investigated in this analysis is ∼2.7 ns.

For clarity, only simulation results from select simulation setups are shown in all the figures, and more examples are provided in the supplementary material.⁵⁸ Particular focus will be on the BLYP–D3–TZV2P simulations as they are the ones with the largest employed basis set and plane wave cutoff and thus should be largely converged with respect to the electronic structure settings for the BLYP–D3 functional approximation. In the discussion that follows, the two sets of BLYP–D3 simulations can uniquely be identified by the choice of Gaussian basis set: DZVP vs. TZV2P. For each simulation setup, the pairs of trajectories referenced below will be denoted by “−1” or “−2” appended to the shorthand name of the simulation settings (e.g., BLYP–D3–TZV2P–1). The temperature of interest will similarly be specified using the NVT target temperature.

In the discussion of proton transfer in recent literature,^31,32 the term “concerted hop” is often said to be important for proton migration, but its exact meaning is often not defined with respect to the time scale that makes a set of multiple hops concerted. A proton hop as defined in this work is the process by which an excess proton is identified as “belonging” to different water molecules in two different time frames of a trajectory. Therefore, the definition of a hop is inherently associated with a time scale τ. One common choice is to use the lower limit of τ for analysis corresponding to the frequency that trajectory snapshots are written, which for our simulations is 0.5 fs. If τ is larger than the period of hopping between two water molecules (i.e., “rattling”), fewer (and possibly much longer) hops are recorded. If τ is too small, too much emphasis can be placed on analyzing rattling events as opposed to actual long-range proton migration. Instead of studying a single τ, the large battery of AIMD simulations in this work can be used to investigate the hopping rate for different processes over a range of τ values. To better understand the relevant time scales for several types of proton hop events, a hop function h_n(τ) is defined to count the number of successful proton hopping events, each of which involves n distinct water molecules in n consecutive time intervals of τ (Figure 1(a)). A more detailed explanation of h_n(τ) is given in the supplementary material. This hopping function is normalized to count the number of successful proton hopping events per 50 ps and its value is independent of the actual path taken by the proton within τ.

The number of successful hopping events, h_n(τ), is plotted in Figure 1(b) for n equal to 1, 2, and 3 and τ spanning 0.5 fs to 0.5 ps. As one should expect, h_n(τ) is sensitive to the time scale used to probe whether a proton transfer event did in fact occur and the length of the water wire involved. Even though proton transfer events along water wires containing four or more water molecules (n > 4) are well-defined, the present AIMD simulations are not of sufficient length to adequately sample those rarer processes. The corresponding curves for four AIMD simulations (the results from the other trajectories are provided in Figure S1 in the supplementary material) display similar features for n = 1 and n = 2 over the entire range of τ sampled, which is encouraging to see some degree of reproducibility across density functional approximations and simulation settings. The situation, however, begins to differ for the case of the longer water-wire hops involving n = 3 water molecules, where now the height and the location of the maximum hop rate peak can be quite different even between runs with similar simulation settings (e.g., BLYP–D3–TZV2P). The apparent discrepancy for the BLYP–D3–TZV2P trajectories is likely a result of insufficient statistics, but the differences shown in these two particular trajectories emphasize their different burst-and-rest behaviors that will be further discussed in relation to other computed observables.

One important question to ask is how these concerted hops along water wires (n = 3) are related to long-range proton migration, as the hop function alone does not convey information about the net displacement (or its direction) over some time interval. In Figure 2, the distance of hydronium oxygen (O^*) relative to its initial position is plotted as a function of simulation time for the same BLYP–D3–TZV2P, BLYP–D3–DZVP–1, and HCTH–2 trajectories. The analysis is general to all other trajectories, and more examples of the burst-and-rest behavior can be found in Figure S2 in the supplementary material. The first thing to notice is that the variance of the distance a proton can travel can be large, such as the case of the BLYP–D3–TZV2P trajectories in which one of the two trajectories exhibits a single burst period and the other does not. To see how concerted hops affect proton migration, the trajectory paths are marked with black empty squares to show the occurrences of all concerted hops with n = 3 and τ = 50 fs as defined in h_n(τ). One can clearly see that even though many of these concerted hop events happen when the proton is migrating, they also appear when the proton is essentially trapped. This says that concerted hops do not always induce proton migration. Furthermore, the frequency of the concerted hops is clearly shown to be quite dependent on the density functional; there are many more concerted hops for the BLYP trajectories (BLYP–2 not shown) than those of BLYP–D3–TZV2P. Similar findings are observed for other choices of τ. Even over the course of 100 ps of multiple AIMD trajectories, these simulations are not yet long enough to observe several periods of burst dynamics separated by rest dynamics. It should be noted that the BLYP functional was used in the recent AIMD study of concerted proton hopping.³¹ 

The burst-and-rest behavior in proton dynamics, as shown in Figure 2, is known from prior literature.^15,31,32 To understand this behavior in more detail, the “burst” and “rest” regions from Figure 2 were separately analyzed to identify possible correlations between these regions and the local H-bond network solvating the excess proton. Recent work by Hassanali et al. provided an interesting analysis using ring statistics^31,59–61 to study the distributions of the number of H-bonded rings and its relationship to the presence of long-lived (lifetime >5 ps) and short-lived traps (lifetime <2 ps). Basically, this analysis counts the number of shortest closed-loop paths, each consisting of the hydronium oxygen and two of its adjacent oxygen neighbors in an undirected graph whose vertices are the oxygen atoms and the edges are the hydrogen bonds. They found in their BLYP simulations, using the DZVP basis set and plane wave cutoff of 300 Ry, that long-lived proton traps tend to be larger at smaller ring count values, whereas the short-lived traps tend to be larger at larger ring count values. The rationale is that fewer rings can stabilize the location of the excess proton (trapping), and vice versa. A similar analysis^62–65 is applied to the “burst” segments that consist of short-lived traps and the “rest” segments that contain long-lived traps. Similar trends are also found in the current BLYP–D3–TZV2P and BLYP–DZVP (see Figures 3(a) and 3(b)) trajectories here, but these trends do not hold here for the BLYP and HCTH simulations (Figures 3(c) and 3(d)). Not only does this suggest that these kinds of analyses can be sensitive to the initial conditions of a trajectory and/or choice of the density functional but it also suggests that a burst period can occur even when the number of water wire pathways from the hydronium is limited.

In Figure 4, the RDF or (g(r)) between hydronium oxygen (O^*) and water hydrogen (H_w) was separately computed for the “rest” and “burst” trajectory segments shown in Figure 2. Noticeable differences in peak positions and intensities in the RDFs are observed from the “burst” and “rest” periods. The water hydrogens in the second solvation shell are further away from the hydronium oxygen for “burst” compared to the “rest” period for BLYP with and without dispersion corrections. For the BLYP functional, the RDF peak for the first solvation shell is shifted to shorter distances during the “burst” period, but no noticeable shift is observed in the RDF for the BLYP–D3 functional. This difference between density functionals may be related to the increased number of n = 3 water-wire hops (19 vs. 5 detected n = 3 events) observed in the “burst” region, suggesting an enhanced lifetime for a compressed water-wire in the BLYP simulation relative to BLYP–D3. Most notable, however, is the density for the RDF peak near 2.0 Å, which is significantly higher for “burst” dynamics. The contribution to this peak comes from water molecules donating a weak H-bond to the hydronium cation, a solvation structure detail known to be highly correlated with successful proton transfer events according to pre-solvation concepts (see Figure 5(a)).^10,28,43 In the discussion that follows, these weak H-bond donating water molecules will simply be referred to as the “fourth water” or the “fourth water neighbor.” Note that the O^*–H_w RDF calculated from the two entire BLYP–D3–TZV2P trajectories that are shown in Figure 2 shows similar differences as the “burst” and “rest” segments. More examples of the O^*–H_w RDFs are provided in Figure S2 in the supplementary material.

Despite its seemingly small contribution in these pair-distribution functions, the role of this fourth water molecule is pivotal in the observed dynamical properties of the excess proton, as emphasized in studies on proton transfer in bulk liquid.²⁸ Details on the excess proton solvation structure involving this fourth water are further investigated here in analysis of joint probability distributions of the distance between the fourth water and hydronium (d) and the angle of approach for the fourth water to the hydronium (θ), as defined in Figure 5(a). The “burst” and “rest” segments of the BLYP–D3–TZV2P–2 trajectory are analyzed separately in Figures 5(b) and 5(c). Two critical features can be observed by comparing these two segments. The first is that during the “burst” period, the fourth water molecule is largely 2.0 Å from the hydronium oxygen (consistent with the RDFs) and populating a cone roughly 10°–20° relative to the plane of the hydronium cation. This angle is not likely to be close to zero.^66,67 Additionally, the overall shape of the two probability distributions shows evidence of contrasting behavior. As the fourth water molecule participates strongly in proton transport during the “burst” period, there is preferential sampling of short distances and angles. During the “rest” period of the trajectory, the water hydrogen loosely stays at 3.0 Å from the hydronium and spans a much larger range of θ angles from 20°–60°. At the larger θ angles relative to the hydronium normal during the “rest” period, there is more competition for the fourth water to donate its H-bond to a nearby water molecule, thus forming a void region (see below) above the hydronium cation. The general shapes of these joint probability distributions for the “burst” and “rest” periods are consistent across all AIMD simulations examined in this work. The strong contrasting relationship in these joint probability distributions for the “burst” and “rest” segments of a trajectory suggests that the distance and angle coordinates employed can be used as good indicators for detecting whether an excess proton is migrating through the H-bond network or trapped as a distorted Eigen complex.⁹ 

The next question to ask is whether the fourth water just happens to be in the right spot during proton migration or if it must first translate and orient itself to enter the “void” cavity residing above the hydronium oxygen (hydrophobic face).^10,20 By Figure 5(b), the fourth water hydrogen is defined to be proximal to the hydronium oxygen if it has a distance d smaller than 2.6 Å and angle θ smaller than 35°. Otherwise, a void exists above the hydronium oxygen. Utilizing these criteria, one can compute the fraction of time that a void is filled ϕ relative to the lifetime of that hydronium before the excess proton hops. In Figure 6(a), a value of ϕ = 0 indicates the void was present during the entire lifetime of the hydronium cation, while ϕ = 1 indicates a void never formed or it was always filled by a fourth water. As learned from the analysis above, it is not surprising to observe that the longest-lasting hydronium centers are those when the void lasts for (almost) the entire lifetime, and when the void is filled, the lifetimes tend to be shorter. As is clear from Figure 6(b), the probability distribution for the fraction of time that the void exists is bimodal, i.e., the void likely either exists for the entire lifetime of the hydronium center or the fourth water is proximal the entire lifetime. This indicates that the fourth water is most likely donating its H-bond to hydronium right from the start or there is competition with a nearby water molecule for this fourth water’s H-bond for most of the lifetime of the hydronium center, thus requiring a fluctuation of the H-bond network to swing down the fourth water’s H-bond towards the hydronium cation facilitating the next hop. It is important to consider that some (or all) of the observed slow behavior is likely related to the known glassy nature of DFT water simulations.

The fraction of time the void is filled ϕ is markedly different when the proton is in the “burst” or “rest” period. Defining the “filled”, or simply “F”, state to be ϕ ∈ [1 − Δϕ, 1] that corresponds to the top bar in the histogram in Figure 6(b), Figure 6(c) reveals that the hydronium is more likely to be in its F state in the “burst” than in the “rest” period. This also suggests that it is more probable for the F state to be retained in the “burst” period after each hop. Evidence for the continuation of the filled state comes from the conditional probability, P(ϕ_n ∈ F|ϕ_n-1 ∈ F), that is defined to be the probability that the current hydronium is in the F state given the previous hydronium was in the F state. The conditional probability being higher in the “burst” period than “rest” period confirms the continuation of the F state is more likely in the “burst” period. This higher likelihood of being continuously in the “filled” state is consistent with the fact that the proton is able to migrate a longer distance in the “burst” period. Figure 6 was prepared using the 300 K BLYP–D3–TZV2P–2 data, but similar trends were seen in the other simulation setups, and a summary is provided in Table S1 in the supplementary material.

One can learn much about the atomistic details of proton transfer by studying how the configurations of atoms in the first two solvation shells around hydronium are correlated to the distance between the fourth water and hydronium. In Figure 7(a), a schematic diagram of the neighboring H-bonded water molecules is shown for the first and second solvation shells. The notation δ(a, b) will be used to indicate the distance between two atoms a and b. For example, the definition of d from Figure 5 would be annotated as δ(4, O^*) defined using the Figure 7(a) schematic.

Figures 7(b) and 7(c) show two examples of how joint probability distributions were used to deduce the correlated motions by neighboring water molecules as the fourth water hydrogen approaches O^*. The black arrows in Figures 7(b) and 7(c) are guides to the eye indicating how the distance along the x-axis is correlated to decreasing d. For example, as d decreases (represented by the direction of the red arrow in Figure 7(a)), Figure 7(b) indicates that the most likely motion for the “1x” oxygen is to approach O^* and this is represented by the direction of the green arrow at “1x” oxygen in Figure 7(a). Similarly, Figure 7(c) implies that the correlated motion for the “H” hydrogen is to move away from O^* as represented by the direction of the blue arrow in Figure 7(a). Following the same procedure, the most likely shifts resulting from approach of the fourth water to O^* can be mapped out. However, it is important to note the directions of arrows drawn in Figure 7 are only most probable, but the shifts are not guaranteed as such because of thermal fluctuations. Second, only a pair correlation between d and a single distance δ is considered at any one time, as convergence of the complete joint probability distribution of d with all distances would require exponentially more data for each dimension to achieve a similar quality of statistics. Therefore, Figure 7(a) addresses the relationships that the red arrow has with each single blue/green arrow. Inferring the relationship between a blue arrow and a green arrow from Figure 7(a) can lead to wrong conclusions. However, one can deduce d is positively (negatively) correlated to the sum of two other distances if both of them are positively (negatively) correlated to d. For example, d is positively correlated to δ(O^*, H) – δ(O^*, 1x), which is commonly used as a proton transfer coordinate, from Figures 7(b) and 7(c). A few other similar joint probability density plots are shown in Figure S3 in supplementary material.

In Figure 7(a), except that of the 4_1x oxygen, all oxygen (green) arrows point inwards toward O^*. This is consistent with the previous observations of a compressed water wire being involved in the transport of an excess proton.^31,68 Indeed, there is an entire solvation sphere of influence contracting onto the O^* as the fourth water approaches O^*. The strength of this contraction becomes weaker as the distance of a water molecule increases away from O^*. Therefore, the effect on the “1x” wing is strongest, then the “1y” wing, and finally the weakest effect is observed on the 1z wing (Figure 7).⁹ It is also observed that the hydrogen (blue) arrows tend to point away from O^* and this is consistent with an increased likelihood for successful proton transfer, which is promoted by approach of the fourth water molecule. Local fluctuations of hydrogen bonding in the solvation shells around a hydronium have been shown to facilitate proton transfer,^8,69 and this effect is likely affected by the presence (or absence) of the fourth water molecule. For the case of the 1z water hydrogens, no arrows were indicated because the distance δ(1z, O^*) was found to be sufficiently far away that only negligible correlations with d were observed for those hydrogens. This is consistent with the finding by Reishcl, Köfinger, and Dellago⁷⁰ that shows the electric field contribution dies after 7 Å. From Figure 4, by observing that the peak at 3.2 Å peak is slightly moved to smaller distances, it hints that the water oxygens are pulled inwards towards O^* more than the first shell hydrogens are pushed away from the O^*.

Additionally, Figure 7(a) indicates that if the 4_1x water is proximal to the soon-to-be hydronium cation, it tends to be slightly repelled as the proton transfer from O^* to the 1x oxygen progresses, i.e., both the 4_1x oxygen and hydrogen atoms move away as d decreases. This is presumably due to the enhanced hydrophobic character of the “1x” water molecule as it accepts handoff of the excess proton (special-pair dance).⁹ This correlation is consistent with previous reports that the coordination number of the 1x water must decrease during proton transfer (pre-solvation theory).^2,4,5,28 Note also that the 4_1x water was previously identified as an important part of the so-called Eigen-Zundel-Eigen (EZE)^2,9,12,13 mechanism for proton transfer, even though the EZE mechanism does not emphasize the role of the fourth water.

As previously reported, the O_w–O_w RDF (Figure 8(a)) computed from DFT approximations as used in these AIMD simulations tends to be over-structured when compared with experiment.^71,72 This over-structuring, although still present, is less severe for HCTH/120 compared to the other density functionals examined here. While the inclusion of exact Hartree-Fock exchange has been shown to be important for modeling aqueous systems (e.g., vibrational properties, electronic structure, solvation of anions),^74,75 the results in this case for the hybrid functional B3LYP are not significantly improved over other generalized gradient approximation density functionals examined here. Some softening of the water-water structure is observed relative to the BLYP and BLYP-D3 results; however, better agreement with the experimental curves amongst the functionals examined here is observed with the HCTH/120 functional at 300 K. In preliminary studies, the M06–L functional⁷⁶ for water-water properties appears to further improve the results of functional approximations examined here, even though the proton-water properties (Figures S4–S10 in the supplementary material) have a larger discrepancy with the experiment as well as the free energy barrier (Figure S15 in the supplementary material) for proton transfer is quite different from the other AIMD setups. Different degrees of improvement were observed when the temperature of the simulations is elevated and the change is most significant for the BLYP–D3–DZVP and HCTH setups.

Comparing the RDFs involving hydronium-water pairs of particles, the elevated temperature does not seem to have a big effect on the O^*–O_w, O^*–H_w, H^*–O_w RDFs (Figures 8(b)–8(d)). Interestingly, there is no density at 2.0 Å for O^*–H_w from the experiment and modeling by Botti, Ricci, and Soper.⁷³ The Empirical Potential Structure Refinement (EPSR) procedure^77–79 that was used to fit the scattering data might have been affected by the use of rigid hydronium and water models. Relative to their flexible counterparts, the rigid models may have provided insufficient flexibility to account for the distorted solvation structures expected of the excess proton associated with rattling events, the special-pair dance,⁹ and approach of the fourth water molecule. It would be particularly interesting to examine the RDFs that result from an EPSR fit to the scattering data using fully flexible and/or reactive molecular models to look for the presence of density near 2.0 Å in the O^*–H_w RDF. Additionally, with the inclusion of a more flexible model, density in the O^*–H_w RDF is expected to decrease (shifting towards smaller separations) to develop a more noticeable minima near 3.7 Å based on comparison of RDFs computed from nonreactive and reactive models, which would bring the experimental RDF into better agreement with simulation below 3.7 Å.¹⁶ Though, to be fair, the solvation structures obtained from these AIMD simulations are too structured relative to the experimental curves, particularly O^*–H_w in both the first (∼3.2 Å) and second (∼5.0 Å) solvation shells. Inclusion of nuclear quantum effects (which is what classical AIMD simulations at elevated temperatures are suggested to mimic) is expected to soften the structure of the solvation shells, which may bring the DFT results into better agreement with experiment. If, in the experimental RDF for O^*–H_w near 2.0 Å, there truly is no density present, then that would suggest that the fourth water molecules contributing to that density are too long-lived in these AIMD simulations. To further improve agreement with experiment and reduce density near 2.0 Å, it would seem that the solvation structure from the AIMD simulations above the hydronium cation would need to be altered in such a way to simultaneously reduce the fraction of time the fourth water spends actually donating its H-bond to hydronium, but still participate in the proton transfer process. More accurate density functional approximations (with nuclear quantum effects) capable of better reproducing the experimental water-water RDFs, such as RPBE-D3,⁸⁰ may facilitate mobility of water molecules solvating the excess proton and contribute to both decreasing the density in the peak near 2.0 Å and increasing the fraction of successful proton transfer events. Larger versions of Figure 8, the O_w–H_w, H_w–H_w, and H^*–H_w RDFs are provided in supplementary material (Figures S4–S10) to facilitate a more careful comparison and study by the reader.

The self-diffusion constant was extracted in the usual manner via the Einstein relation as 1/6 of the slope of mean-squared displacement (MSD) as a function of time in its linear regime. Please refer to Figures S11–S14 in supplementary material for the corresponding plots of MSD for each AIMD simulation. The estimated self-diffusion constants for both water and excess proton are shown in Figure 9. At 300 K, the relatively small values for the BLYP and BLYP–D3 water diffusion coefficients are in agreement with previous reports.^27,29 The water diffusion coefficients determined from the HCTH simulations are about twice smaller than previous simulations computed using the Car-Parrinello method for shorter times.^26,83,84 Increasing the temperature from 300 K to 330 K (actual NVE temperatures reported in Table I) roughly increases the self-diffusion constant of water by a factor of 3 (Figure 9(a)), whereas the proton self-diffusion constants are less sensitive to the temperature change (note the change in scale of the y-axis in the Figure 9(b) plots), with the exception of the HCTH setup in which the diffusion constant increased by a factor of about 3.7 (Figure 9(b)). Consistent with previous discussion, the proton diffusion constants for the BLYP–D3–TZV2P trajectories vary widely between the runs with not having an identifiable “burst” period during the first simulation. The diffusion constant for the second simulation, which did contain a “burst” period, is the largest determined from the current set of AIMD simulations at 300 K. This result along with the (expected) increased variability of the proton diffusion constants relative to the water diffusion constants is a strong indicator that multiple long-time simulations (more than that computed in this work) are necessary to statistically converge dynamical properties to at least a couple of significant figures.

Considering the slow water dynamics and the fact that no quantum nuclear effects were taken into account, the overall proton diffusion is surprisingly fast from these AIMD calculations. Indeed, by using the obtained averages from each pair of simulations, the following ratios of the excess proton self-diffusion constant to that of water were determined: 31 (BLYP–D3–TZV2P), 70 (BLYP), 55 (BLYP–D3–DZVP), and 11 (HCTH), whereas the ratio from experiments is 4 at 300 K. Increasing the temperature to 330 K in the AIMD simulations improves agreement of these ratios with experiment at 300 K but does not completely resolve the issues with the over-structured and slowly diffusing DFT water liquid. The ratios from AIMD at 330 K were determined to be 14 (BLYP–D3–TZV2P), 31 (BLYP), 15 (BLYP– D3–DZVP), and 11 (HCTH).

There are two main components to the proton transport mechanism: vehicular and hopping. The vehicular component is known^12,85,86 to contribute a smaller percentage to the overall proton diffusion, thus a factor of 2 or 5 times slower translational water diffusion and corresponding decrease of the vehicular contribution would not alone significantly decrease the magnitude of the proton diffusion coefficients. However, a result from this decreased water diffusion may be that intact water wires remain stable for long times relative to the length of these short AIMD simulations providing prime conditions for the protons to rapidly translocate along H-bonds. Moreover, it is well known these DFT approximations generally give lower barriers for proton transfer between water molecules compared to higher-level quantum calculations.^87,88 Conversely, the relatively slow water dynamics could increase the lifetimes of voids and local proton traps preventing successful long-range proton transport. The case of the BLYP–D3–TZV2P simulations at 300 K is one example of these two possible scenarios with no “burst” period in the first trajectory and the second trajectory possessing a “burst” period with 5 water-wire hops (n = 3) detected and the proton transported nearly 22 Å in 25 ps (Figure 2(a)). The “burst” period of the BLYP–1 trajectory (Figure 2(b)) is an example where a relatively large number of water–wire hops was detected (19), but the overall distance traveled by the proton was only 6–8 Å after the initial successful proton hop effectively stalling proton transport until a solvent fluctuation initiated the onset of a “rest” period. Note that such a solvent fluctuation could have easily initiated a “rest” period while the proton was near its original location at the start of the “burst” period. The relatively short time scales accessible in these AIMD simulations can make it challenging to fully address the statistical properties of proton transport and its coupling to solvent fluctuations in the surrounding H-bond network. The typical DFT approximations used in the AIMD simulations also give relatively poor properties for water and proton transfer barriers, but these two errors appear to actually cancel to some degree to give seemingly accurate results for proton diffusion relative to experiment. However, if one divides the AIMD proton diffusion result by the AIMD water self-diffusion constant and compares that ratio to experiment, the agreement is quite poor. Nevertheless, the AIMD simulation data can still be used to derive computationally efficient models with better bulk properties that can extend the results of these important AIMD simulations and explore the mechanistic details of proton transport in multi-nanosecond trajectories.^16,84

The term “concerted hop” can be ambiguous because its meaning is often implied in the context, but not from a definition. An important question investigated in this work is whether clearly defined concerted hops always result in long-range proton migration. From a sizeable suite of AIMD simulations with various density functional approximations and temperatures, it was determined that even though concerted hops, as defined with respect to the time scale τ, are important for long-range proton migration, concerted hops can also be observed even when the proton is in a “rest” state. Furthermore, the frequency of concerted hops involving three participating water molecules over time scales τ in the range of 0.5 fs to 0.5 ps was found to be dependent on the density functional, basis set, and initial conditions. One of the BLYP–D3–TZV2P trajectories in this study had the lowest frequency of concerted hops while the identity of the simulation with the highest frequency of concerted hops was found to vary as a function of τ. The computed AIMD trajectories (most longer than 100 ps) were further analyzed with the identification and partitioning into “burst” and “rest” segments. From this analysis, a fourth water molecule, which donates a weak H-bond to hydronium, tends to have a much stronger presence in “burst” periods than in “rest” periods as determined through analysis of joint probability density distributions and the presence of a key signature peak near 2.0 Å in the O^*–H_w RDF. Absence of density for this RDF peak was correlated with long-lived traps for the proton. Conversely, this fourth water hydrogen being proximal to the hydronium oxygen promotes proton migration consistent with the concepts of pre-solvation theory.

An important question examined in this work is whether the fourth water is always in a prepared state to donate its H-bond to hydronium or if it must first translate and orient itself from outside of the void, possibly facilitated by solvent fluctuations, to donate its H-bond to the hydronium oxygen. Statistics from the large ensemble of AIMD simulations (totaling ∼2.7 ns) indicate that in about 80% of the proton hopping events analyzed, the fourth water is either already there from the start donating a H-bond or it is not close enough to donate its H-bond for the entire lifetime of the new hydronium center. In the latter scenario, the fourth water is unlikely to move into the hydrophobic void before the next successful proton hop. This leads to the hypothesis that the observed burst/rest activities of the excess hops are caused by a series of hops onto different hydronium centers that are or are not already pre-solvated by the fourth water neighbor due to statistical fluctuations of the water H-bond network. In effect, the fourth water neighbors of the hydronium cation and dynamic candidate water for being protonated (i.e., the 1x water in Figure 7) are involved in a game of tug-o-war influenced by their respective solvation shells. Similar to previous work, water-wire compression was also observed here and found to be modulated by at least the first two solvation shells. Finally, analysis of the changes in structural and dynamic properties for two different temperatures determined that the hydronium-water RDFs were not as sensitive to temperature as the water O_w–O_w RDF. Notably, the experimental RDF by Botti et al. determined using the EPSR procedure did not yield any density near 2.0 Å in the O^*–H_w RDF. While it is not clear if this discrepancy is caused by the simulations or by the analysis of the scattering data, it is evident that this particular signature of the O^*–H_w RDF with contributions from the fourth waters needs to be understood completely in the future because of its importance to proton dynamics. It is conjectured that a combination of a more flexible model in the EPSR fit and use of a more accurate electronic structure method (e.g., density functional) will meet in the middle and converge onto similar radial distribution functions with density in the O^*–H_w RDF near 2.0 Å. Considering the critical role of these fourth water neighbors, special attention should be given when developing force fields for studying proton transfer in aqueous environments. Preliminary results along this route with the development of a new reactive MS-EVB model do in fact indicate that the excess proton self-diffusion is increased with enhanced density near 2.0 Å in the O^*–H_w RDF.

This research was supported in part by the National Science Foundation (NSF), through Grant No. CHE-1214087. Portions of this work were supported by (U.S.) Department of Energy (DOE) under Contract No. DE-AC02-06CH11357. S.T. acknowledges the Croucher Foundation for a postdoctoral research fellowship. The computational resources in this work were provided in part by a grant of computer time from the U.S. Department of Defense (DOD) High Performance Computing Modernization Program at the Engineer Research and Development Center (ERDC) and Navy DOD Supercomputing Resource Centers and in part by the University of Chicago Research Computing Center (RCC).