Proton mobility in aqueous systems: combining ab initio accuracy with millisecond timescales

Abstract

We present a multiscale simulation of proton transport in liquid water, combining ab initio molecular dynamics simulations with force-field ensemble averaging and kinetic Monte-Carlo simulations. This unique Ansatz allows for ab initio accuracy incorporating the femtosecond dielectric relaxation dynamics of the aqueous hydrogen bonding network, and bridges the time-scale gap towards the explicit simulation of millisecond diffusion dynamics.

---

1 Introduction

Proton transfer in water plays an important role in many chemical and biological processes.^1–3 Especially with regard to Proton Exchange Membrane (PEM) fuel cells, the high conductivity of water is a desirable property and the reason for the efficiency of Nafion-based membranes.^4–8 These membranes enclose water channels, in which proton conduction is possible. In this context, a lot of research has been conducted to further clarify the proton conduction mechanisms both experimentally and theoretically.

A model picture of the underlying mechanisms which explain the anomalously high proton conduction in water was actually first proposed by Grotthuss more than 200 years ago.⁹ It comprises two main elements: a hopping step, which is later followed by reorientation. The hopping step can actually happen as a series of proton transfers along a chain of water molecules, or rather “water wire”,^10,11 which leads to the propagation of the excess charge along the chain without any diffusive motion of the molecules themselves, and a simultaneous reversal of the chain's dipole moment. Consequently, the subsequent rotation of the water molecules is necessary in order to enable further proton transfer along the same “water wire”.

In recent years, this idealized model picture has been refined further. In particular, the water reorientation by rotation is believed to be an oversimplification as a water molecule in bulk water is hydrogen-bonded to four neighbor molecules, which hinders free rotation. Instead, the breaking of a hydrogen bond, which is necessary to reduce the hydrogen-bond coordination number of an accepting H₂O from four to three, is considered the rate-determining step.^12,13

With the continuously increasing computing capabilities, ab initio calculations of proton transfer in water have become feasible. It is now possible to investigate proton transfer on the atomistic level.^13–19 However, while offering a high accuracy, these methods are still computationally expensive and give only access to timescales of tens to hundreds of picoseconds. Especially with regard to the simulation of macroscopic properties of proton exchange membranes, they often do not allow statistically relevant sampling of the excess proton movement in a reasonable amount of time. Similarly, the adequate incorporation of co-solvents, geometric confinement, effects, or the influence of functional interfaces is computationally challenging.

In this work, we present a coupled Molecular Dynamics/Lattice Monte Carlo (cMD/LMC) approach, which aims at simulating large-scale excess charge diffusion in aqueous systems. This combined method has already been applied successfully to PEM materials^20,21 and solid acids.²² Based on structural dynamics calculated from molecular dynamics (MD), it builds a time-dependent excess charge transition matrix for the Lattice Monte Carlo (LMC) scheme. The mapping between the topology of the system and the resulting excess charge transfer rate is determined with statistical analysis from ab initio MD.

On the basis of a neutral classical MD trajectory of water, we show that our cMD/LMC scheme is able to determine dynamical quantities of the excess charge in best agreement with ab initio methods over large time scales.

2 Method

For the propagation of an excess proton through a liquid water system, we have developed a hybrid scheme which combines aspects from molecular dynamics (MD), and lattice/kinetic Monte Carlo methods.^23–25 Protons can jump between lattice sites with transfer rates which are extracted on-the-fly from an underlying MD simulation (see also ref. 20–22). Each lattice site represents an oxygen atom of the considered system and can be in one of two states: occupied or unoccupied. The trajectory-based evolution of the pairwise oxygen distances in the MD results in a change of the proton transition rate in the LMC scheme. This scheme has been shown to yield a realistic picture of condensed-phase systems.

In order to take the characteristic hydrogen-bond network of H₃O⁺ in liquid water into account, we construct the LMC topology such that each oxygen site is connected to its closest three neighbors.

2.1 Lattice Monte Carlo scheme

The LMC system consists of a discrete range of states and is described by a master equation with time-dependent transition rates:where p_i(t) is the probability to find the system in state i, and ω_ij is the transition rate from state i to state j.

The scheme first translates the transfer rate from a lattice site i into a “lifetime” of this very site. In this context, the lifetime is defined as the duration τ, after which the excess charge at lattice site i is transferred to an adjacent site. In our specific algorithm, this lifetime is calculated by drawing a random number x in the interval (0, 1] and solving the equation

for τ.

Once this specific lifetime has passed, the algorithm determines the target site for the excess charge by drawing another random number x from the interval and choosing the second oxygen site with index j₀ such that

For a more detailed mathematical explanation, see Prados et al.²³

2.2 Determination of the excess charge diffusion

The challenge in determining the excess charge motion from an atomic trajectory is the fact that there is no “labeled” excess charge particle whose position can be tracked. Instead, similar to the Grotthuss picture, a multitude of protons are involved during a (long-range) excess charge transfer.

Within the ab initio molecular dynamics (AIMD) simulations, the excess charge diffusion is therefore determined by means of a collective variable specifically designed for this task, with an effective coordinate given by:

R_O and R_H in eqn (4) stand for the positions of the oxygen atoms and hydrogen atoms, respectively. R^excess has the property that it is invariant under translations of neutral H₂O molecules while being sensitive to the motion of single (excess) protons and H₃O⁺ complexes. This effective variable can easily be adjusted for more complex chemical systems (as has been demonstrated for R-PH₂O₃/R-PH₃O₃²⁶).

In contrast to MD, the cMD/LMC scheme represents the excess charge position explicitly as an overprotonated molecular state that is assigned to a water molecule at each time step. It is therefore trivial to follow the excess charge over time and determine its diffusion.

Unlike MD, the proton motion in the cMD/LMC scheme is not continuous: when a water molecule transfers its excess charge to another nearby molecule, the excess charge position changes instantaneously from the position of the donor molecule to the position of the acceptor molecule. As a consequence, the mean squared displacement function is discontinuous and changes stepwise by the full oxygen–oxygen distance.

2.3 Excess charge transfer rate determination

In order to obtain a realistic transfer rate matrix between individual lattice sites for the cMD/LMC scheme, the instantaneous atomic geometry around a protonated water molecule must be mapped onto the probability of an excess charge transfer for a given constellation of the H₃O⁺'s surroundings. We determine this mapping by statistically analysing the excess charge transfer events within the ab initio MD trajectory and determining the transfer rate as a function of the oxygen–oxygen distance between H₃O⁺ and the surrounding H₂O molecules by means of the following algorithm for each time step t:

(1) Determine the distance distribution ρ^OO(d^OO,t) of potential H⁺ acceptor water molecules around the H₃O⁺.

(2) Lookup the distance distribution ρ^jump(d^OO,t) of actual H⁺ acceptor molecules, i.e. consider only those waters to which a proton jump has actually occurred in the next MD step at t + Δt.

(3) The proton jump rate at time t is determined via

The final jump rate is obtained by averaging over all time steps

ω^jump(d^OO) = 〈ω^jump(d^OO,t)〉_t (6)

We fit this jump rate to a function of the formThe resulting fit parameters, from our ab initio molecular dynamics simulation are listed in the ESI.†

3 Results

3.1 Analysis of ab initio molecular dynamics of protonated water

3.1.1 Trajectory setup. For this work, ab initio trajectories of protonated H₂O and a classical trajectory of neutral H₂O were simulated. We used the software CP2K^27,28 for the ab initio simulations, and Gromacs^29–31 for the classical MD simulation.

While the ab initio trajectories are used both as a reference system and for analysis of the dynamical properties of H₃O⁺, the classical trajectory provides the underlying topology for calculating the excess charge jump rates of the cMD/LMC scheme.

AIMD trajectories of 100 water molecules and an excess proton in a periodic box of length 14.405 × 14.405 × 14.405 ų were used to determine the proton jump rates and analyze the radial distribution function (RDF) between hydronium and neutral water. In order to investigate the temperature dependence of the observables, three trajectories at 300 K, 350 K and 400 K were run. We used the BLYP functional^32,33 with a DZVP-MOLOPT-SR-GTH basis and the GTH pseudopotential.³⁴ The DFT-D3 dispersion correction was used.³⁵ The molecular dynamics simulation was run with a timestep of 0.5 fs using the Nosé–Hoover thermostat.³⁶

At the force field level, an MD trajectory consisting of 216 water molecules in a periodic box of length 18.626 × 18.626 × 18.626 ų was run on a simulation time scale of 10 ns at 350 K. For the classical MD, we used the TIP4P water model with the Amber force field.³⁷

3.1.2 Ab initio excess charge diffusion. We determine the mean squared displacement (Fig. 1) from R(t) (see Section 2.2). After about 5 to 10 ps, the mean squared displacement (MSD) reaches an almost linear regime. The remaining fluctuations are due to the limited statistics (only one excess proton is present in the water box). From a linear fit we get a diffusion coefficient D = (1.48 ± 0.01) Ų ps⁻¹. We use this value as a reference to judge the results of the cMD/LMC approach.

Fig. 1 Mean squared displacement of the excess charge within the ab initio trajectory. The slope of the linear fit is used to determine the diffusion coefficient.

3.1.3 Excess charge transfer rate. We determine the excess charge transfer rate according to Section 2.3 from the available AIMD trajectories at 300 K, 350 K and 400 K. The resulting transfer rates as functions of the donor–acceptor distance are shown in the ESI.† We observe a noisy development of the jump rates at the borders of the distance axis. This is due to the fact that only few O–O distances in these regions are sampled within the AIMD trajectories. In the middle region, on the other hand, the better sampling leads to a more steady course. We take the varying sampling into account by weighting the jump rate results with the number of observed hydronium–H₂O pairs N^OO. All three jump rates yield very similar results, which is in agreement with previous studies.²⁰ In the following sections, the trajectory at 400 K will be used.

3.2 Dielectrical relaxation following an elementary proton hopping step

The presence of an excess proton in water leads to considerable electrostatic forces, which influence the structure of water molecules in the vicinity of the H₃O⁺ ion. This effect is incorporated automatically by AIMD, which is known to yield a precise description of both structure and femto/picosecond dynamics around a solvated charge,^15,16,38 but is missing within the cMD/LMC scheme, as the excess charge state within the LMC is not returned as feedback to the MD simulation. As a result, the acceptor–donor distances in the LMC scheme are overestimated, which leads to a diminished proton hopping rate and therefore an underestimation of the overall excess charge diffusion.

We propose a scheme which adjusts the trajectory data to take into account this dielectric response of the water upon protonation changes. Specifically, we analyze how the presence of an excess charge in the AIMD influences the RDF of the first solvation shell. This change is subsequently used to construct a conversion function that provides a mapping between the pairwise oxygen distances of an uncharged water system and those of a protonated system within the first solvation shell.

3.2.1 Construction of a conversion function. The microsolvation of an H₃O⁺ cation differs in two ways from that of a neutral H₂O molecule: first, the three donor hydrogen bonds are shorter, and second, the acceptor bridge is elongated. Since our force-field MD only contains H₂O molecules, we have to consider the different solvation structure within the Monte Carlo part, in order to yield more realistic excess charge transfer rates.

To address this issue, we construct a conversion function which maps the distance distribution of an uncharged water and its first solvation shell to the distance distribution of a hydronium ion and its first solvation shell. This allows us to take the structural effect of a charged particle within the LMC scheme into account by effectively rescaling the distances between an occupied oxygen site and its neighbors.

For the construction of the conversion function, we use the integrated RDF, which yields the number of particles within a radius r. It is defined as

with the particle density ρ and the radial distribution function g(x).

We also define n_hydronium(r) as the integrated RDF of the desired distribution (i.e. the hydronium–H₂O distribution), and n_neutral(r) as the integrated RDF of the original distribution in the neutral H₂O system.

The conversion function c(d) is then defined as:

c(d) = n_hydronium⁻¹(n_neutral(d)) (9)

using the inverse function of n_hydronium. The lower picture of Fig. 3 shows the resulting conversion.

In our implementation, we calculate the conversion function at discrete points, and use a linear interpolation scheme to provide a continuous mapping.

Note that the conversion function effectively shortens the donor hydrogen bonds by around 0.25 Å in the range between 2.5 Å and 3 Å. This agrees well with the distance of the peaks of the H₃O⁺–H₂O versus the H₂O–H₂O (see Fig. 2).

Fig. 2 Comparison of the radial distribution functions between neutral water molecules, and between a hydronium ion and water. The dotted blue line with a peak at 2.5 Å shows the radial distribution function of a hydronium ion and water. The dashed green line and the dashdotted red line show the radial distribution functions of water from our ab initio and classical MD simulations.

Fig. 3 Upper picture: Cumulative distribution functions of the oxygen–oxygen distribution for neutral water, and for neutral water and a hydronium ion. The arrows indicate how, given an unrescaled distance of a neutral system d^OO_neutral, one arrives at the corresponding rescaled distance d^OO_hydronium. Lower picture: Resulting conversion function, which maps a distance d^OO_neutral from the oxygen distance distribution of uncharged water molecules onto the distance d^OO_hydronium from a distance distribution between an H₃O⁺ and a neutral water molecule.

3.2.2 Modeling relaxation time within the LMC scheme. Upon the jump of a proton from one water to another, the hydrogen bond network responds in the way described above (Section 3.2.1). However, this response is not instantaneous, but has a characteristic relaxation time. In order to simulate this relaxation of the water molecules in the vicinity of the newly emerged H₃O⁺ ion, we rescale the O–O distances when an excess proton has just jumped (at time t₀) linearly in time from the unrescaled (d^u) to the rescaled (d^r) distances.with t₀ < t < t₀ + t_relax.

The choice of the dielectric relaxation time parameter t_relax determines the time span, in which the linear rescaling takes place. It turns out that an adequate range for t_relax is found between 2 and 10 ps.

3.3 Proton dynamics in the cMD/LMC approach

In the following, we investigate how the cMD/LMC scheme reproduces dynamical quantities such as MSD and the excess charge autocorrelation. For this, we compare the dynamical properties of the LMC scheme with the reference AIMD calculation and experimental results.

3.3.1 Diffusion coefficients. Table 1 lists the resulting diffusion coefficients of the LMC scheme with different parameters, the reference AIMD trajectory, and experimental data. Ab initio and experimental results differ by about a factor of two. This is a good agreement, as routinely factors between 0.1 and 10 are reported.^40,41

Table 1 Comparison of the diffusion coefficients of the cMD/LMC scheme and a reference AIMD simulation

Method	D/Å^2 ps^−1
LMC (w/rescaling, trelax = 0 ps)	10.5(3)
LMC (w/rescaling, trelax = 4 ps)	1.69(6)
LMC (w/rescaling, trelax = 6 ps)	1.72(5)
LMC (w/o rescaling, (trelax → ∞))	1.01(3)
Reference AIMD (this work)	1.48(1)
Experimental	0.865^39

---

Next, we compare the results of the cMD/LMC scheme, and ab initio results. Fig. 4 shows the variation of the excess charge diffusion coefficient with respect to the relaxation time parameter. We see that the diffusion coefficient is in best agreement with the reference AIMD for values t_relax between 4 and around 12 ps. Values of t_relax below 2 ps, on the other hand, lead to unrealistically high diffusion coefficients which differ by an order of magnitude. Using the oxygen distances of the neutral water system yields a lower diffusion coefficient (around 1 Ų ps⁻¹) than our reference AIMD. The reason for this behavior can be seen from Fig. 2: the neutral oxygen–oxygen distances lie around a distance of 1.7 Å, which corresponds to a proton jump rate of less than one jump per ps. It should be noted that the lowest possible diffusion rate in our scheme occurs, if the excess charge does not jump at all. In this case, only the oxygen diffusion contributes with D = 0.6 Ų ps⁻¹.

Fig. 4 Diffusion coefficients of the cMD/LMC scheme against the relaxation time parameter t_relax(x). The dashdotted line shows the diffusion coefficient of the cMD/LMC scheme with unrescaled distance (representing the limiting case for t_relax → ∞). The dashed line represents the diffusion coefficient of the ab initio simulation at 1.48 Ų ps⁻¹ as a reference.

While the concept of the relaxation time parameter t_relax is physically motivated, the range of 4 to 12 ps obtained by comparison of the final diffusion coefficient to the reference AIMD simulation may appear a quasi-empirical choice. However, this range can also be derived from an analytical consideration: assuming an activation energy of 2.6 kcal mol⁻¹ for the cleavage and reformation of two hydrogen bonds,⁴² and using a wavenumber of 170 cm⁻¹ for the hydrogen bond stretching frequency,⁴³ application of the Arrhenius equation yields an average time of around five picoseconds before a hydrogen bond breaks. This quantitative agreement validates the choice of t_relax = 5 ps as a universal value (for liquid water).

3.3.2 Excess state lifetime. A complementary parameter for the characterization of the dynamics of the excess proton is the lifetime of an (excess) protonation state of a water molecule. In order to compare our LMC scheme against our reference AIMD simulation, we first define a state function i⁺(t), which yields the oxygen index of the current H₃O⁺ ion at time t. We then define the excess charge autocorrelation function η as

η(t) = 〈δ_{i+(t0),i+(t0+t)}〉_t0 (11)

where δ is the Kronecker delta. The excess charge autocorrelation function can be interpreted as the probability that an oxygen atom that was threefold protonated at time t₀ is still threefold protonated at time t₀ + t.

Fig. 5 shows the results for the reference AIMD and LMC the scheme for a range of values for t_relax. Qualitatively, the results are similar to the diffusion coefficients: the best agreement with the reference ab initio trajectory is found for values of t_relax between 2 and 8 ps. The LMC scheme with unrescaled distances shows the lowest decline of η, caused by the larger oxygen–oxygen distances and the resulting decreased excess charge transfer rates. The LMC scheme with instantaneous relaxation, on the other hand, shows the fastest decline.

Fig. 5 Autocorrelation of the H₃O⁺ state for different values of t_relax. The upper triangular line shows the autocorrelation of the cMD/LMC scheme with unrescaled distances as the limiting case for t_relax → ∞. The dotted line shows the autocorrelation for t_relax = 0. The solid line shows the autocorrelation of the ab initio simulation as a reference.

4 Conclusion

We have investigated proton dynamics in liquid water on extended time scales using a multiscale scheme which combines ab initio and classical molecular dynamics with kinetic Monte-Carlo. By virtue of deriving all relevant energetic/kinetic parameters needed for the KMC part from ab initio calculations, we are able to reach an accuracy comparable to quantum chemical methods on mesoscopic time scales.

The diffusion constants and the excess state lifetimes obtained from our multiscale Ansatz agree quantitatively with reference AIMD simulations. Our studies show that the dielectrical relaxation of the immediate surroundings of the hydronium ion needs to be explicitly incorporated into the scheme for a realistic description of the excess proton diffusion process. The typical relaxation time in our scheme agrees well with hydrogen bond lifetimes, which we see as a confirmation of the importance of hydrogen bond cleavage as the rate determining step for excess proton transfer.

Our approach is easily applicable to more complex multi-component systems such as Nafion, water channels in proteins, and other nano structured molecular systems.

Conflicts of interest

There are no conflicts to declare.

References

1. M. Rini, B.-Z. Magnes, E. Pines and E. T. J. Nibbering, Science, 2003, 301, 349–352.
2. O. F. Mohammed, D. Pines, J. Dreyer, E. Pines and E. T. J. Nibbering, Science, 2005, 310, 83–86.
3. S. Meiboom, J. Chem. Phys., 1961, 34, 375–388.
4. S. Dokmaisrijan and E. Spohr, J. Mol. Liq., 2006, 129, 92.
5. M. K. Petersen and G. A. Voth, J. Phys. Chem. B, 2006, 110, 18594–18600.
6. D. Seeliger, C. Hartnig and E. Spohr, Electrochim. Acta, 2005, 50, 4234–4240.
7. S. Feng and G. A. Voth, J. Phys. Chem. B, 2011, 115, 5903–5912.
8. M. Cappadonia, Solid State Ionics, 1995, 77, 65–69.
9. C. J. T. De Grotthuss, Mémoire sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique, 1805.
10. G. Bekçioğlu, C. Allolio and D. Sebastiani, J. Phys. Chem. B, 2015, 119, 4053–4060.
11. G. Bekçioğlu, C. Allolio, M. Ekimova, E. T. J. Nibbering and D. Sebastiani, Phys. Chem. Chem. Phys., 2014, 16, 13047–13051.
12. N. Agmon, Chem. Phys. Lett., 1995, 244, 456–462.
13. D. Marx, ChemPhysChem, 2006, 7, 1848–1870.
14. L. Vilčiauskas, M. E. Tuckerman, J. P. Melchior, G. Bester and K.-D. Kreuer, Solid State Ionics, 2013, 252, 34–39.
15. M. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Chem. Phys., 1995, 103, 150–161.
16. M. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Phys. Chem., 1995, 99, 5749–5752.
17. M. E. Tuckerman, D. Marx and M. Parrinello, Nature, 2002, 417, 925–929.
18. D. Marx, M. E. Tuckerman, J. Hutter and M. Parrinello, Nature, 1999, 397, 601–604.
19. P. Intharathep, A. Tongraar and K. Sagarik, J. Comput. Chem., 2006, 27, 1723–1732.
20. G. Kabbe, C. Wehmeyer and D. Sebastiani, J. Chem. Theory Comput., 2014, 10, 4221–4228.
21. G. Kabbe, C. Dreßler and D. Sebastiani, J. Phys. Chem. C, 2016, 120, 19905–19912.
22. C. Dreßler, G. Kabbe and D. Sebastiani, J. Phys. Chem. C, 2016, 120, 19913–19922.
23. A. Prados, J. J. Brey and B. Sánchez-Rey, J. Stat. Phys., 1997, 89, 709–734.
24. Y. Cao, H. Li and L. Petzold, J. Chem. Phys., 2004, 121, 4059–4067.
25. D. T. Gillespie, J. Phys. Chem., 1977, 81, 2340–2361.
26. G. A. Ludueña, T. D. Kühne and D. Sebastiani, Chem. Mater., 2011, 23, 1424–1429.
27. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
28. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
29. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
30. D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
31. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
32. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
33. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
34. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
35. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
36. D. J. Evans and B. L. Holian, J. Chem. Phys., 1985, 83, 4069–4074.
37. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
38. C. Allolio and D. Sebastiani, Phys. Chem. Chem. Phys., 2011, 13, 16395–16403.
39. J. H. Simpson and H. Y. Carr, Phys. Rev., 1958, 111, 1201–1202.
40. J. C. Grossman, E. Schwegler, E. W. Draeger, F. Gygi and G. Galli, J. Chem. Phys., 2004, 120, 300–311.
41. E. Schwegler, J. C. Grossman, F. Gygi and G. Galli, J. Chem. Phys., 2004, 121, 5400–5409.
42. S. Cukierman, Biophys. J., 2000, 78, 1825–1834.
43. W. B. Bosma, L. E. Fried and S. Mukamel, J. Chem. Phys., 1993, 98, 4413–4421.

---

Footnote

† Electronic supplementary information (ESI) available: Proton jump rate functions determined from AIMD, illustration of the linear distance rescaling. See DOI: 10.1039/c7cp05632j

---