Accurate diffusion coefficients of the excess proton and hydroxide in water via extensive ab initio simulations with different schemes

[ This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in The Journal of Chemical Physics, 157, 024504 (2022) and may be found at: https://doi.org/10.1063/5.0093958 ] Despite its simple molecular formula, obtaining an accurate in silico description of water is far from straightforward. Many of its very peculiar properties are quite elusive, and in particular obtaining good estimations of the diffusion coefficients of the solvated proton and hydroxide at a reasonable computational cost has been an unsolved challenge until now. Here I present extensive results of several unusually long ab initio MD simulations employing different combinations of the Born--Oppenheimer and 2$^\mathrm{nd}$-generation Car--Parrinello MD propagation methods with different ensembles (NVE and NVT) and thermostats, which show that these methods together with the RPBE-D3 functional provide a very accurate estimation of the diffusion coefficients of the solvated H$_3$O$^+$ and OH$^-$ ions, together with an extremely accurate description of several properties of neutral water (such as the structure of the liquid and its diffusion and shear viscosity coefficients). In addition, I show that the estimations of $D_\mathrm{H_3O^+}$ and $D_\mathrm{OH^-}$ depend dramatically on the simulation length, being necessary to reach timescales in the order of hundreds of picoseconds to obtain reliable results.


I. INTRODUCTION
One of the most important and intriguing phenomena occurring in water is the diffusion of its two autodissociation products, namely the solvated proton and hydroxide ions. As it is reflected in classical textbooks on electrolytes, 1 the first measurements of their abnormally high diffusion coefficients were seen with fascination, and it was clear since then that the explanation of this fact would be intimately linked to the nature of water itself. 2 While it is now accepted that the diffusion of these species occurs via the Grotthuss mechanism, 3 which involves consecutive proton transfers between neighboring molecules, even after many decades of intense investigation the particular details of such structural migration mechanism and the nature of the solvated proton and hydroxide are still a matter of debate. [4][5][6][7][8][9][10][11][12] As seen in those studies, interpreting the experimental results is not straightforward, and it is necessary to employ computer simulations to obtain a complete picture of the diffusion process of these species. However, despite the vast computational work already done, 3,[13][14][15] an accurate and affordable method capable of reproducing the experimental diffusion coefficients of the H + (aq) and OH − (aq) was yet to be found. Indeed, water is one of the most difficult substances to accurately model in silico, and even finding a reliable approach for simulating pure water is being a challenge in itself. In particular, a) Electronic mail: daniel.munozsan@upm.es within the ab initio molecular dynamics (AIMD) community there is a long-standing and ongoing quest for an ideal exchange-correlation (XC) functional that is capable of reproducing all of its most important properties. 16 The performance of a given XC functional for liquid water is usually assessed by evaluating static properties related to the structure of the liquid, like its radial distribution functions g OO (r), g OH (r) and g HH (r), and the equilibrium density ρ at a given temperature and pressure, with dynamical properties such as the diffusion coefficient D H2O in second place. Ideally, obtaining a proper computational description of the H + (aq) and OH − (aq) implies finding an approach that reproduces their experimental diffusion coefficients in addition to reproducing the diffusion coefficient of neat water. Thus, it is worth to start the search of such optimal approach for H + (aq) and OH − (aq) by reviewing the performance of the different XC functionals for reproducing the experimental diffusivity of H 2 O. For a long time, ab initio simulation of water was mostly done employing bare GGA XC functionals, mainly BLYP 17,18 and PBE 19 (for a detailed review on these and other functionals for simulations of water, including comprehensive tables summarizing results that will be commented in the following, see Ref. 16). A nagging problem of PBE -and of BLYP to a lesser extentknown since long ago 20 is that it overstructures water, and several simulations using indistinctly light or heavy water have established that its diffusion coefficients D around 300 K are in the order of 0.2-0.6 ×10 −9 m 2 s −1 for BLYP and 0.1-0.3 ×10 −9 m 2 s −1 for PBE, 16 roughly one order of magnitude smaller than the experimental values of both H 2 O and D 2 O (respectively 2.39 × 10 −9 m 2 s −1 and 1.96 × 10 −9 m 2 s −1 at 300 K and 1 bar 21 ). This issue has been usually remedied by increasing the simulation temperature, as suggested since the very discovery of the problem. 20 The development of different dispersion correction schemes represented a very good way of improving the GGA results at a negligible computational cost. While adding dispersion corrections in different forms to the BLYP functional greatly improves its description of water and remediates the overstructuring problem 16 (though the reported D for heavy water is at most about 1.7 × 10 −9 m 2 s −1 at 308 K or even higher T , thus still somewhat below the experimental value), the situation is more intricate for PBE. The DRSLL dispersion correction scheme 22 does soften the structure of PBE water [23][24][25] and improves the D D2O values towards 1.7-2.1×10 −9 m 2 s −1 at 301-304 K, 23,24 while TS only improves modestly the radial distribution functions without increasing much the low diffusivity (D D2O = 0.44 × 10 −9 m 2 s −1 at 300 K, 26 ) and no improvement in either is obtained with DCACP 27 corrections 28 or with the popular D2 28 and D3 corrections. 25 Actually, the fact that D2/D3 corrections are contraindicated for PBE simulations of water was already demonstrated at the cluster level in the seminal paper where Grimme et al. introduced the D3 correction scheme, 29 as seen in Fig. 8 therein where PBE-D2/D3 are shown to perform much worse than bare PBE for the WATER27 benchmark set.
Despite the PBE functional is still one of the most used for simulations of water, probably because its robustness for a wide range of systems makes it a good choice for simulations of aqueous interfaces and heterogeneous systems, in recent years its two modified versions revPBE 30 and RPBE 31 (often including dispersion corrections and/or a fraction of Hartree-Fock exchange) are becoming new favorites for simulating water. Interestingly, revPBE and RPBE are the only GGA functionals that do not produce a severe overstructuring of the liquid. 16 The bare revPBE functional yields D D2O = 2.1 × 10 −9 m 2 s −1 at T = 323 K, 28 which is in very good agreement with the experiment when taking into account the finite-size effects of the simulation. Other simulations in heavy water show that adding DCACP 28 or D3 25 corrections worsens the diffusivity by reducing it with respect to both the bare revPBE and experimental values, while D2 28 or DRSLL 23 worsen it by its increase. On the other hand, more recent simulations 32 with light water show that revPBE-D3 yields a very good (finite sizecorrected) value of D H2O = 2.22 × 10 −9 m 2 s −1 at 300 K, while an even better estimation (2.29 × 10 −9 m 2 s −1 ) and a remarkably exact description of the structural and spectral features of water are obtained when using revPBE0-D3 (which adds a 25% of Hartree-Fock exchange over revPBE-D3) and including nuclear quantum effects (NQEs) via path integral MD. About the importance of including NQEs in molecular simulations of water, it is worth noting that years ago these were thought to slightly weaken the H-bonds (and thus increase the diffusivity) on the basis of experimental evidence such as isotope effects in light vs. heavy water. This motivated approximations such as slightly elevating the temperature in order to compensate for the missing NQEs in simulations with classical nuclei. However, in recent years it has emerged a more complex picture where NQEs originate competing effects which both weaken and strengthen the H-bonds, 33,34 depending on the precise features of the aqueous system and on the property of interest whether these cancel out or else one of them dominates. At that respect, the results obtained with revPBE-D3 and revPBE0-D3 with and without NQEs 32 show that adding NQEs to revPBE-D3 reduces the diffusivity (from 2.22×10 −9 m 2 s −1 to 1.60×10 −9 m 2 s −1 ), but when added to revPBE0-D3 they are largely compensated by the fraction of exact (HF) exchange which has the opposite effect, in such a way that bare revPBE-D3 (which includes neither ) still provides a very good description of water. This improved understanding of the impact of NQEs discourages the use of crude approximations such as elevating the temperature, which may be actually counterproductive. 35,36 On the other hand, the effect of NQEs on simulations differs depending on the underlying method. For instance, simulations with neural network potentials trained with the SCAN meta-GGA functional 37 showed that the D H2O = 1.1 × 10 −9 m 2 s −1 with classical nuclei slightly decreased to 1.08 × 10 −9 m 2 s −1 when including NQEs at 300 K (but more noticeably at higher temperatures); but simulations with the meta-GGA functional B97M-rV 38,39 showed the opposite effect 35 (increase from D H2O ∼ 2.7 × 10 −9 m 2 s −1 to ∼ 2.9 × 10 −9 m 2 s −1 at 300 K, values inferred from the plot in Fig. 3 of Ref. 35). I note in passing that meta-GGA functionals do not perform necessarily better than regular GGA functionals for water, as shown by the previous examples.
Finally, the RPBE functional is nowadays mostly used together with D3 corrections. RPBE-D3 performs very well for different aqueous systems, from pure bulk water [40][41][42][43] and aqueous solutions [44][45][46] to interfacial 47 and nanoconfined water, [48][49][50] including the liquid/vapor equilibrium 51,52 and the liquid at supercritical conditions. 51,53-55 The diffusivity of RPBE-D3 water was first estimated via neural network potentials, 56 which reported D H2O in the range 3.0-3.2×10 −9 m 2 s −1 for systems of different sizes at 300 K and ρ = 1 g/cm 3 after applying the proper finite-size correction (see Fig. S4 in the SI Appendix of Ref. 56). More recently, it was also estimated directly via AIMD with a quite thorough protocol involving averaging over several tens of NVE simulations, 53 yielding D H2O = 2.0 × 10 −9 m 2 s −1 . This disagreement between the direct AIMD estimation and the neural network one is very likely due to the underlying training set employed for the neural network potentials.
After having observed the notorious difficulties for nailing down the experimental D of neutral water via AIMD, it should come as no surprise that obtaining accurate estimations of D H3O + and D OH − has been a challenging task. For these, on top of the aforementioned problem of finding a good functional, there is the problem of the extremely long simulation times required to properly converge the diffusion coefficients in comparison to neutral water: obviously, a simulation with N neutral water molecules yields a statistics for D H2O that is N times that of a simulation of the excess proton/hydroxide, where only one of these species per unit cell is usually employed.
Pioneering AIMD simulations 57 with heavy water showed that the bare BLYP functional reproduces qualitatively the experimental D D3O + /D OD − ratio, though their absolute values are below the experimental ones by roughly a factor of 2. On the other hand, other functionals like HCTH and PW91 produce qualitatively wrong results when comparing D D3O + to D OD − . 57 Interestingly, those three functionals yield qualitatively similar values of D D3O + and also of D D2O , but differ greatly in the estimation of D OD − which varies more than order one of magnitude. The underlying reason is that these different functionals produce fundamentally different solvation structures of the hydroxide ion, which in turn results in qualitatively different structural migration mechanisms. 13,57 Very recent simulations 12 using BLYP-D3 and two different empirically-corrected versions of it, analyzing two replica systems in each case, yielded different values of D H3O + in a wide range (2.9-9.8 × 10 −9 m 2 s −1 ). However, there are no data for the OH − (aq) with the same methods.
The best AIMD estimation so far of D for the solvated proton and hydroxide was recently obtained using the PBE0-TS functional and employing heavy water at 330 K, 7 obtaining D D3O + = 8.3 × 10 −9 m 2 s −1 and D OD − = 3.7 × 10 −9 m 2 s −1 , in fair agreement with the respective experimental values in heavy water at 298 K of 6.7-6.9×10 −9 m 2 s −1 and 3.1-3.2×10 −9 m 2 s −1 . The observed underlying migration mechanisms of the hydronium and hydroxide basically agree with those proposed years back, 13 though pointing out now a prominent role of concerted proton transfer events in the case of H + migration, an issue which is however debated. 12 Despite such overall good results, it must be noted that PBE0-TS somewhat overestimates the experimental D for the hydronium and hydroxide while it underestimates that of neutral water even when using a simulation temperature of 330 K 26 (1.46 × 10 −9 m 2 s −1 in the D 2 O simulation at 330 K vs. the experimental 1.96 × 10 −9 m 2 s −1 of heavy water at 300 K). Such simulation temperature is used in Refs. 26 and 7 in order to approximately compensate for the missing NQEs, 26 which implies that their computed D values must be compared to the experimental ones at 300 K. However, it must be stressed that this approximation may actually produce the opposite effect compared to the actual NQEs, 32,36 as said above.
Regarding the importance of NQEs for modeling the diffusion processes of H 3 O + and OH − in water, it must be realized first that these are governed by the hydrogen-bond fluctuations of the water molecules in their second and first solvation shells, respectively, and not by the energetic barriers for proton transfer between the H 3 O + /OH − and their solvating waters. 13 According to this picture, NQEs must affect H 3 O + /OH − diffusion mostly through the indirect effect they have in the overall dynamics of water. Indeed, seminal path integral AIMD simulations 58,59 showed that including NQEs did not change these overall diffusion mechanisms 13,57 but only changed some peculiar features of both defects. In the case of the proton, NQEs facilitate the interconversion of the idealized Zundel and Eigen complexes (resp. 58 which as said is not the limiting step for the structural diffusion. For the hydroxide, NQEs favor proton transfer within the transient (H 2 O) 3 OH − complex for a wider range of geometries compared to the classical case where a more stringent geometrical arrangement is required; 59 this could slightly favor productive proton transfer (i.e. structural diffusion) without qualitatively changing the mechanism. 13 Looking at the isotope effects experimentally observed for the diffusion coefficients (0.82 for D 2 O/H 2 O, 0.72 for D 3 O + /H 3 O + and 0.59 for OD − /OH − , all at 301.2 K; 21,60 note that isotope effects in dynamical properties are not entirely due to NQEs but also partly due to trivial effects because of the greater masses of the substituted atoms that appear also in the classical limit 33 ), these observations may explain why these do not differ much between neutral water and hydronium and differ a bit more in the hydroxide case. Taking all this into account, it is reasonable to expect that NQEs would not change much the relative values D H3O + /D OH − , D H3O + /D H2O and D OH − /D H2O , while predicting their impact on the different D values in absolute terms seems mere speculation at this point.
All the information laid out above makes it clear that there is still quite room left for quantitatively improving the computational estimations of the diffusion coefficients of H 3 O + and OH − , hopefully without requiring expensive approaches such as adding HF exchange or explicitly including NQEs. Motivated by the good results in the literature of RPBE-D3 for many properties of water in very different situations, and in particular for its diffusion coefficient at ambient conditions, I performed extensive AIMD simulations to compute the diffusion coefficients of the excess proton and hydroxide ions in bulk water at the RPBE-D3 level. More precisely, the two main goals reached in this paper are: (i) providing reliable estimations of D H3O + and D OH − by performing long and accurate Born-Oppenheimer MD simulations, and (ii) showing that it is possible to obtain equally reliable estimations of D H3O + and D OH − at a much smaller computational cost exploiting the 2 nd -generation Car-Parrinello MD scheme.

II. METHODS
All the simulations were performed with CP2K. 61 Unless otherwise stated, the model system for neutral water consisted in 128 H 2 O molecules in a cubic cell with L = 15.663 Å (resulting in a water density of 996.556 kg/m 3 , matching the experimental value at T = 300 K and p = 1 bar 62 ), in which a proton was added or removed to obtain the water systems with a H 3 O + or OH − , respectively. The electronic structure of the system was treated within the spin-restricted formulation of Kohn-Sham theory with the RPBE functional 31 plus D3 corrections 29 which included only the two-body terms and had zero damping. I employed GTH pseudopotentials [63][64][65] together with a TZV2P basis set, 66 using a 500 Ry density cutoff for the expansion of the auxiliary planewave basis set and applying NN50 smoothing for both the charge density and its derivatives. As shown in the literature, 67 NN50 smoothing reduces the spurious numerical effects related to the evaluation of the exchangecorrelation energy on a finite grid, and tests done in the present study showed that it improves the accuracy of the forces in a way equivalent to a cutoff increase of at least 100 Ry. In addition, tests for determining the accuracy of the TZV2P basis set in comparison with other basis sets available for CP2K and with a planewave basis in the near-infinite limit (see Sec. I of the Supporting Information (SI)) showed that TZV2P is probably the best choice in terms of performance vs. accuracy for liquid water, being almost as accurate as QZV3P and better than the popular DZVP-MOLOPT-SR set. 68 Two kinds of ab initio MD simulations were carried out: Born-Oppenheimer (BO) 69 and 2 nd -generation Car-Parrinello (2genCP) 70,71 MD. In the BOMD runs, I set a threshold for the SCF convergence of 2 × 10 −7 and a timestep of 0.4 fs; the latter allows an excellent accuracy for integrating the equations of motion while also facilitates the SCF convergence w.r.t. larger timesteps as e.g. 0.5 fs (compensating to some extent the cost of performing more MD steps for the same simulation time). The BOMD simulations were carried out either in the NVE ensemble or in the NVT ensemble; in the latter case with a Nosé-Hoover chains thermostat 72 with a chain length of 3, an order of 9 for the Yoshida integrator and a time constant of 100 fs. On the other hand, the 2genCPMD runs were carried out either using its original formulation which implies a Langevin-like thermostatting, or using other thermostat schemes (Nosé-Hoover chains or CSVR 73 ). In all cases, I used a 0.4 fs timestep for consistency with the BOMD settings. The full details will be unfolded in the Results section.
In all cases, the diffusion coefficients were determined as usual via the Einstein's relation MSD(t) = 6Dt where the mean square displacement (MSD) was computed by following the position of the O atom of the molecule of interest (in the case of neutral water, averaging the MSD over all the water molecules in the system). The plots of every MSD(t) function and the time ranges used for computing the corresponding slopes 6D for all the diffusion coefficients in this work are disclosed in Sec. III of the SI.

III. RESULTS AND DISCUSSION
A. Determination of finite-size effects and estimation of the viscosity As known since long ago, 74,75 diffusion coefficients in liquid systems estimated via MD employing periodic boundary conditions are affected by finite-size effects. In the present case of a cubic simulation cell, the 'exact' coefficient D 0 in an infinite system and the diffusion coefficient actually obtained in the simulation D PBC are related by the expression D 0 = D PBC + kBT ξ 6πηL , being k B the Boltzmann constant, T the temperature, the numerical factor ξ ≈ 2.837297, η the shear viscosity of the solvent, and L the length of the simulation cell.
A convenient way of determining the finite-size correction term kBT ξ 6πηL consists in computing D PBC for different system sizes. The subsequent linear representation of D PBC vs. 1/L readily provides D 0 as the y-intercept and the correction term as the slope, from which in turn the shear viscosity η of the liquid can be estimated. It must be noted that a priori there is no guarantee that this value of the viscosity will be the same as the experimental one, and therefore using the latter to directly estimate the correction term kBT ξ 6πηL may imply further errors. On the other hand, MD estimations of η do not show any significant system-size dependence. 75 To this aim, three BOMD runs of neutral water systems were carried out employing respectively 64, 128 and 256 water molecules in cells of sizes 12.431 Å, 15.663 Å, and 19.734 Å. The three were carried out in the NVT ensemble with a Nosé-Hoover chains thermostat with the settings detailed in the Methods section. All were equilibrated during 10 ps, followed by production runs of 800 ps, 400 ps and 200 ps respectively for the 64, 128 and 256 H 2 O cases.
As shown in Fig. 1, the dependence of D PBC on 1/L is linear as expected, providing an estimation of D 0 = 2.26 × 10 −9 m 2 s −1 which is in excellent agreement with the experimental D exp = 2.39×10 −9 m 2 s −1 at T = 300 K. 21 Moreover, the estimation of the shear viscosity that is readily obtained is η sim = 7.3 × 10 −4 Pa · s, which is in remarkably good agreement with the experimental value η exp = 8.9 × 10 −4 Pa · s. 77 The resulting additive correction term for the simulations containing 128 water molecules is 0.546 × 10 −9 m 2 s −1 at T = 300 K. In the following, all the reported diffusion coefficients D H2O , D H3O + and D OH − include this correction term kBT ξ 6πηL (where T is the average temperature of each simulation).
The previous excellent estimations of the dynamical properties D and η of H 2 O are accompanied by an ac- In the previous section I have estimated a timedependent property such as D via MD simulations in the NVT ensemble, as it is also usually done in the literature. However, it must be kept in mind that the very role of the employed thermostat is 'tweaking' the dynamics of the system in order to maintain the canonical sampling, in such a way that a priori there is not any guarantee that the resulting dynamics will provide the proper time-dependent properties. The most accurate way of estimating such observables is to generate a set of initial conditions (i.e. positions and momenta) representative of a canonical ensemble, and with them run many (usually in the order of several tens) microcanonical simulations that will indeed provide the proper dynamics of the system, as done in e.g. Ref. 53 (see Fig. 2 therein for a very clarifying scheme of such simulation protocol). Averaging the results of that set of simulations will yield the proper estimation of the desired observable for the canonical ensemble. However, the amount and length of the simulations required to properly converge this ensemble average make such procedure cumbersome when done with AIMD, which is why time-dependent properties are often computed via a single NVT run or sometimes even a single NVE run.
It is frequently assumed that the diffusion coefficients obtained via NVT simulations are acceptable (except in pathological cases such as when using very strong stochastic thermostats, as I will discuss later), and indeed in some cases these are remarkably accurate, as it is the case of the NVT results for D H2O shown in the previous section. However, in order to firmly demonstrate that an estimation via NVT is also justified for the particular cases of acidic and basic water, and that the NVT and NVE results are consistent, I performed 1 NVT and 5 NVE runs for each system (including also the neutral water system for the sake of consistency). While 5 NVE runs are of course too few to obtain a proper ensemble average, they are enough to make a semi-quantitative comparison between the NVE and NVT approaches and discern whether they are in agreement.
In order to obtain the initial conditions for the different NVE runs, after an initial equilibration of 20 ps I carried out a NVT run with a massive Nosé-Hoover chains thermostat (with the same settings described in Sec. II for the regular (global) Nosé-Hoover chains thermostat), from which 5 sets of positions and momenta where collected at intervals of 5 ps for each system. The use of a massive Nosé-Hoover chains thermostat in this particular case guarantees that the sampled sets of positions and momenta are decoupled from each other even after only 5 ps of time separation. These initial conditions were used to launch NVE runs of 200 ps length for neutral water and 400 ps for acidic and basic water. On the other hand, after an initial equilibration of 10 ps I carried out a NVT run with a global Nosé-Hoover chains thermostat of 400 ps length for each system (note that for the neutral water case, such simulation was already performed during the determination of the finite size effects).
The obtained diffusion coefficients are listed in Table I. The NVE results are dispersed over a certain range that is roughly centered around the NVT results, in such a way that there exists a certain correlation between the average temperature in each simulation and the observed diffusion coefficient D, with a dependence that in the case of H 3 O + and OH − is quite pronounced as visualized by the plot of D vs. T in Fig. 3. This temperature dependence is similar to that shown by the experimental diffusion coefficients D H3O + and D OH − , which is greater than in the case of D H2O . However, despite this general trend, some NVE simulations with similar T may yield different D values, as it occurs for the NVE runs with T ∼ 290 K for both H 3 O + and OH − . Assuming the best-case scenario that different NVE simulations with the same T must yield the same D in the limit of infinite simulation time, 80 these results indicate that 400 ps is quite short to properly converge D H3O + or D OH − via microcanonical sampling. This observation clearly warns against the practice of computing D based on single NVE simulations, as it is sometimes done in the literature.
Overall, even though it is evident that several more NVE simulations would be needed for an accurate determination of the diffusion coefficients via their ensembleaveraging, it seems clear that the NVT-predicted coefficients are consistent with the NVE ones. This conclusion is of course system-and thermostat-dependent, and so far it is only completely valid for the precise thermostat settings used in this case.
Regarding the comparison to the experimental results, the present estimation of D H2O is excellent as already said in Section III A, while the estimations for the hydronium and hydroxide ions are extremely good despite being somewhat below the experimental values. However, I note that the ratio D exp H3O + /D exp OH − = 1.8 at 298 K is equal to the ratio D sim H3O + /D sim OH − provided by the NVT results at 300 K, while the ratio D H3O + /D H2O is 3.9 (experiment) vs. 3.8 (simulation), and the ratio D OH − /D H2O is 2.2 (experiment) vs. 2.1 (simulation). Therefore, the current simulation strategy perfectly reproduces the relative diffusive behavior of all the species H 2 O, H 3 O + and OH − in bulk water at ambient conditions. C. 2 nd -generation CP MD: setting the stage I have just shown that it is possible to obtain excellent estimations of the diffusion coefficients of H 2 O, H 3 O + (aq) and OH − (aq) via Born-Oppenheimer MD in the NVT ensemble, though the simulation times involved are rather high considering the high computational cost of a BOMD simulation. It is then worthwhile to try alternative approaches for performing NVT simulations with Born-Oppenheimer accuracy at a smaller computational cost. This is precisely the case of the 2 nd -generation  Table I Table I, lines between points are a guide to the eye.
Car-Parrinello propagation method, 70,71 where the full SCF procedure during each MD step is substituted by a prediction-correction scheme. At each MD step, this involves the extrapolation of the wavefunction from the previous steps' wavefunctions, followed by one or two correction steps that are equivalent to the usual SCF steps in a regular BOMD. After testing different settings for the stability and performance of the propagation, I chose to employ an order of 2 for the ASPC extrapolation 82 and a fixed number of 2 corrector steps per MD step; I note that the convergence of the wavefunction is not enforced in any further way during the full simulation length. This was enough to obtain a ∼ 2.5× speedup w.r.t. the BOMD simulations. In its original formulation, the 2genCP method allows to perform the canonical sampling of a system by propagating it with a Langevin-type equation: being Ξ I a random noise term defined by Ξ I (0)Ξ I (t) = 6(γ L + γ D )M I k B T δ(t). The reason behind such approach was that regularly integrating the Newtonian equations of motion with the forces obtained via the predictor-corrector scheme F PC led to a dissipative dynamics (though this effect is system-dependent, being less important in systems like water). 70 Thus, by introducing the friction and white noise terms −γ L M IṘI and Ξ I where the latter is defined not by γ L but by (γ L + γ D ), a Langevin dynamics is established where the white noise term is 'enhanced' in order to compensate for the intrinsic dissipative dynamics of the system. In practice, a convenient way of setting up such simulations is choosing a given γ L and then adjusting the value of γ D on-the-fly during the equilibration of the system until obtaining the desired average temperature. This method is therefore an excellent alternative to BOMD when it comes to performing long, expensive MD simulations requiring ab initio accuracy, as it is the case of the determination of diffusion coefficients. However, computing dynamical properties in cases where the underlying dynamics is stochastic -as it is the case of Langevin dynamics-is a delicate issue. While it has been shown that the 2genCP method provides good estima-tions of both static and dynamical properties of systems like water 83 when using a small γ L (∼ 10 −5 fs −1 ), it must be recognized that small values of γ L necessarily imply longer simulation times until the canonical sampling is completely established, while on the other hand higher values of γ L eventually destroy the true dynamics of the system, thus spoiling the estimation of time-dependent properties.
It is therefore desirable to perform a systematic study to find which parameters provide the best sampling of the systems of interest. To that aim, for each of the three studied systems I set up four 2genCP simulations using different values of the γ L parameter (specifically, 0, 10 −5 , 10 −4 , and 10 −3 fs −1 ) starting with γ D = 0 in all cases and equilibrating the systems during 10 ps. Surprisingly, after letting all the simulations to evolve over up to 400 ps of production time, it was observed that the dynamics of all the studied systems was not dissipative at all, even though γ D was in all cases kept as 0 during the whole simulation. Instead, inspection of the running average of the temperature as a function of the simulation time (see plots in Sec. II of the SI) showed that the temperature in the simulations with γ L = 0 tended to stabilize around some value not far from 300 K (similarly to what happened in the BO-NVE simulations), while the simulations with non-zero γ L values behaved as regular NVT runs with different thermostat strengths. On the other hand, the fluctuations of the instantaneous temperature w.r.t. the average value (quantified as the variance σ 2 T , see Sec. II in the SI) in the γ L = 0 cases were similar to those observed in the BO-NVE simulations, while in the γ L > 0 cases they approached the ideal value of σ 2 T for a proper NVT simulation 84 as γ L increased. This behavior was consistently observed for all the systems -neutral, acidic and basic water-. In addition, static properties of water like g OO (r) and g OH (r) were identical to those obtained with BOMD (Fig. 2), while the D values were comparable to the BOMD ones for H 2 O, H 3 O + and OH − , with only small differences that are easily explained by the different values of γ L as will be discussed in the next section.
Thus, apparently the 2genCP method produces in this case a propagation regime that is identical to the BO one without even the need to set a non-zero γ D . The following analyses of the forces and energies during the simulations show that it is indeed the case. Sample structures of the three studied systems were taken every 0.1 ps from the simulations with γ L = 0, and the atomic forces and Kohn-Sham energy of the system were recalculated by completely quenching the wavefunction to the Born-Oppenheimer surface. As shown in Fig. 4, the distribution of the force modules in these structures is practically identical in both the 2genCP and BO cases, while the distribution of the module of the force deviation confirms that the differences between the CP and BO forces are extremely small, practically always well below 5 × 10 −4 a.u. It is interesting to realize that the extent of these force deviations is much smaller than the deviations associated to changing the basis set from the most exact one available in CP2K to the infinite limit within the same level of theory, as seen by comparing Fig. 4 here to Fig. S1 in the SI. On the other hand, the deviation of the Kohn-Sham energy during the 2genCP simulations w.r.t. the BO reference remains always in the order of ∼ 10 −5 Hartrees (Fig. 5), with a remarkably good energy conservation (Fig. 6).
It is then clearly demonstrated that the dynamics provided by the 2genCP propagation scheme at the employed level of theory is indistinguishable from the Born-Oppenheimer regime even with γ D = 0, as the differences between the 2genCP and BO atomic forces are negligible in practice. From this we can conclude that the previous runs with γ L = 0 are actually NVE simulations (note that Eq. 1 reduces to Newton's 2 nd law when γ L = γ D = 0), while those with γ L > 0 are NVT runs with different thermostat strengths as it has been confirmed by analyzing the behavior of the temperature (Sec. II in the SI).
This remarkable behavior of the 2genCP dynamics deserves closing this section with some comments on its relationship to the employed settings. The crucial detail that determines the ability of the scheme for keeping the system so close to the BO surface is the number of corrector steps applied, as stated in the original paper. 70 While therein it is not explicitly investigated whether the dissipative behavior depends on the number of corrector steps, it seems obvious that it must vanish at some point  4. Histogram of the force modules obtained for 4000 snapshots sampled during the 2genCP run of neutral water with Langevin dynamics with γL = γD = 0 (red), in comparison with the corresponding ones obtained at the Born-Oppenheimer surface (black, practically coincident with the red histogram) and with the module of the difference between the 2genCP and BO forces (blue). The inset shows a magnification of the latter. The histograms obtained for the acidic and basic water systems are indistinguishable from these. as the system is kept closer and closer to the BO surface, which in this case occurs with only 2 corrector steps. I note that here the choice of 2 corrector steps was not due to intentionally looking for that effect, but because of its superior accuracy/efficiency balance compared to using 1 corrector step: adding the second step increases the computational cost only a factor of 1.3, but in conjunction with the ASPC extrapolation guarantees an excellent stability of the dynamics. Ensuring such stability was necessary in this work, given the long simulation times pursued, the need for producing a very accurate dynamics in order to properly compute time-dependent properties, and to guarantee that the fast reorganization of the electronic degrees of freedom during the H + /OH − diffusion would remain controlled.

D. Diffusion coefficients via 2 nd -generation CP MD with different thermostats
After observing that in the studied systems the dynamics provided by the 2genCP method is not dissipative, it is straightforward to realize that it is then possible to use it in conjunction with different thermostats, since now the necessity of employing a Langevin equation for the propagation scheme has vanished. Hence, for each of the studied systems I set up three more simulations keeping the previous settings for the 2genCP propagation except the original Langevin thermostat, which is substituted by either a Nosé-Hoover chains thermostat (with the same settings used for the BO-NVT runs), or a CSVR thermostat (with two different values of the time constant, τ = 100 fs and τ = 1000 fs). While Nosé-Hoover chains allow a correct estimation of D as seen above, CSVR 73 can be considered the global version of Langevin (which is local) and therefore it affects the dynamics to a much smaller extent, allowing a better estimation of time-dependent properties. 85 All these runs were equilibrated for 10 ps, followed by 400 ps of production. Together with the previous 2genCP simulations within the original Langevin dynamics scheme, this makes a total of 7 simulations using 2genCP for each system, the results of which I analyze in the following.
The obtained D values are summarized in Table I and visualized in Fig. 3. Overall, the agreement with the diffusion coefficients estimated via BOMD is excellent, with small differences associated to the different thermostats.
The Nosé-Hoover chains thermostat produces the same results for the BO and 2genCP regimes (within the statistical errors), which are also almost equal to those provided by CSVR with τ = 100 fs. The canonical sampling is very well established in these cases, judging by the T fluctuations (Figs. S2-S4 in the SI). This is also the case of CSVR with τ = 1000 fs, which produces the correct temperature fluctuations even though the average temperature slightly deviates from 300 K, and its D values perfectly agree with the previous ones when considering its marked dependence on T . On the other hand, Langevin dynamics with γ L = 10 −5 or 10 −4 fs −1 has more marked problems to produce the required average temperature, being off by about 5 K in at least one of the systems, and with the temperature fluctuations being below the expected value always for γ L = 10 −5 fs −1 and occasionally for γ L = 10 −4 fs −1 . In contrast, Langevin dynamics with γ L = 10 −3 produces both the correct T and σ 2 T in all systems. It is interesting to realize that the D values in the γ L = 10 −5 fs −1 case do nicely agree with those obtained via Nosé-Hoover and CSVR thermostats (in addition to all the BO simulations) when considering the D(T ) dependence, while in the γ L = 10 −3 and 10 −4 fs −1 cases the D H3O + slightly overestimates that of the other approaches. It is then quite possible that the latter cases present the well-known problem of the underlying stochastic dynamics affecting the estimation of time-dependent properties, in this case only observed for the (local) Langevin thermostat.
Finally, the simulations with γ L = 0, which are equivalent to the BO-NVE ones, provide estimations of D which are perfectly in line with those, once more reflecting the correlation between the observed D and T .

E. Analyzing the accuracy of D depending on the simulation length
While it is obviously desirable to achieve the longest simulation times possible, reaching scales of hundreds of picoseconds at the ab initio level as I have done here is a daunting task. As shown above, methods such as 2genCP may alleviate this cost, but even with this approach performing these simulations requires quite some computational effort. The situation would be much worse in the case that one would request an even higher level of theory, such as going beyond the pure GGA approximation to include a fraction of HF exchange or considering nuclear quantum effects. That would surely limit the affordable simulation time to less than 100 ps, probably at most 50 ps, which obviously would very negatively impact on the accuracy of the results. In order to understand such impact, I computed the diffusion coefficients of H 3 O + (aq) and OH − (aq) over different parts of the BOMD simulations with lengths of 50 ps (Fig. 7) and 100 ps (Fig. 8).
As clearly seen in Fig. 7, 50 ps is an extremely short simulation length to reliably compute the diffusion co-efficient of the solvated proton or hydroxide, with the possible results being dispersed over a huge range. Even looking only at the results within the NVT ensemble, which show a smaller dispersion than the NVE ones, it is possible to obtain estimations that are not only quantitatively but also qualitatively wrong (i.e. obtaining D H3O + < D OH − ). Increasing the simulation length to 100 ps slightly improves the situation (Fig. 8), and the estimation of D H3O + via NVT simulations would be now better converged, while that of D OH − could still be rather off. Again, the NVE results show a dispersion that makes impossible to accurately estimate the diffusion coefficients from single NVE simulations, especially if the great dependence of D on T is ignored (I note in passing that the usual procedure of starting a NVE simulation from a state equilibrated via NVT does not at all guarantee that the average simulation temperature during the NVE run will be the same as during the previous equilibration run). This clearly shows the necessity of averaging over at least several tens of NVE simulations  in order to obtain reliable results.

IV. CONCLUSIONS
The results above exposed demonstrate that the combination of the RPBE functional with D3 dispersion corrections provides a very accurate estimation of the diffusion coefficients of the excess proton and hydroxide ions in water, yielding a D H3O + /D OH − ratio that is equal to the experimental one; in addition to an excellent description of different properties of neutral water such as its diffusion coefficient, its shear viscosity, and its liquid structure, in almost perfect agreement with the experimental values. On the other hand, based on the existing results for the similar revPBE-D3 functional, 32 it can be expected that improving these RPBE-D3 results would require adding both NQEs and a fraction of HF exchange, as including only NQEs would produce a (spurious) reduction of D H2O and probably also of D H3O + and D OH − (since their diffusion mechanisms would be now limited by the reduced fluctuations of the H-bonding network in such less mobile water), pushing all three farther from the experimental values. On the other hand, it is difficult to predict how the interplay of HF exchange and NQEs on RPBE-D3 would affect the absolute vales of D H3O + and D OH − . In any case -and even if partly due to such error compensation-, the fact that RPBE-D3 performs so well for all the three species H 2 O, H 3 O + and OH − at an affordable computational cost compared to much more expensive approaches, and without resorting to other approximations usually taken in the field (e.g. increasing the simulation temperature), makes it the best choice so far for simulating them.
Regarding the simulation strategy, the 2 nd -generation Car-Parrinello MD method yields excellent results that are indistinguishable from those obtained with regular Born-Oppenheimer MD at a smaller computational cost. Actually, 2genCP provides in this case a propagation that is so close to the Born-Oppenheimer regime that it is not necessary to employ the Langevin dynamics scheme in which it was originally formulated, being possible to use instead other regular but global thermostats such as Nosé-Hoover chains or CSVR, both of which improve the accuracy of the computed time-dependent properties over Langevin dynamics.
The statistical accuracy of the results presented here is strongly backed-up by the fact that the many simulations performed, employing different propagation regimes and thermostat settings (which guarantees that the different simulations are completely decoupled), provide essentially the same results when taking into account the statistical uncertainties and especially the strong dependence of the diffusion coefficients on the simulation temperature, with only minor differences associated to the different thermostat schemes employed. Moreover, I show the huge importance of reaching simulation times in the order of hundreds of picoseconds when computing the diffusion coefficients of the solvated proton and hydroxide. At that respect, it must be recognized that the times of 400 ps per simulation achieved in this work are much longer than those seen in most of the ab initio MD studies in the literature, where it is usual to find estimations of the diffusion coefficients of these species with simulation times in the order of 50-100 ps. However, as I have just shown, these would be indeed too short to obtain reliable values of D H3O + and D OH − , being rather possible to obtain qualitatively wrong results. Moreover, special care must be taken when computing diffusion coefficients from NVE simulations, since it is necessary to average over at least tens of them in order to obtain reliable estimations, and in particular the great dependence of D on T must be acknowledged.
Finally, I close this paper with some outlook regarding the underlying structural diffusion mechanism of the H 3 O + (aq) and OH − (aq) ions. An initial analysis of the obtained trajectories supports the currently standing picture of both the preferred solvation structures of both ions and their migration mechanisms, i.e. the 'traditional Grotthuss mechanism' 3 in the case of H 3 O + (aq) and the 'dynamical hypercoordination mechanism' 3,13 for OH − (aq). However, deeper analyses of the vast amount of generated data are underway. Such analyses, which are out of the scope of this paper and will be unfolded in following publications, are aimed at understanding more involved aspects of the diffusion mechanism of the excess proton and hydroxide, such as the precise determination of energetic barriers for proton transfer, the occurrence and timescale of concerted proton transfer events, and the landscape of the different solvation structures of both species.

SUPPORTING INFORMATION
The SI document includes a detailed account of the tests for validating the employed basis set, as well as detailed plots showing the evolution of the temperature and its fluctuations during the different simulations, and the MSD computed in all the simulations. In order to check the validity of the basis set, I employed two test systems consisting of 32 water molecules in a cubic unit cell with a = 9.87Å in which a proton was added or removed to obtain respectively an acidic and basic bulk water system. In order to generate an ensemble of representative structures for these systems, two respective AIMD runs were carried out in the NVT ensemble within the 2genCP framework, employing a 0.4 fs timestep and a Nosé-Hoover chains thermostat at 300 K (thus equal to the 2genCP runs with Nosé-Hoover chains described in the main text). The electronic structure of the system was treated as described in the main text for all the systems (i.e. using GTH pseudopotentials together with TZV2P basis sets, a 500 Ry cutoff and using the RPBE-D3 functional). The systems were equilibrated during 10 ps, followed by production runs of 100 ps. Sample structures were taken every 1.0 ps, obtaining 101 structures for each system.  With these sets of structures, the atomic forces were recomputed with CP2K by completely quenching the wavefunction to the Born-Oppenheimer surface using the following basis sets: TZVP, TZV2P, QZV3P, TZVP-MOLOPT, TZV2P-MOLOPT, TZV2PX-MOLOPT, and DZVP-MOLOPT-SR. For this, I kept intact the rest of the electronic structure settings except the D3 dispersion correction, which was removed in order to facilitate the comparison between the different basis sets. In addition, in order to obtain a valid reference for the 'exact' forces within the employed level of theory, I used the planewave code Abinit. This code allows to use the very same GTH pseudopotentials as CP2K, which together with a PW basis set with a cutoff of 900 Rydberg provides in practice the limit of an infinite basis set.
After computing all the individual forces, the appropriateness of the Gaussian basis sets was checked by computing the module of the deviation of the forces with the Gaussian basis w.r.t. the forces computed with the planewave basis, i.e. |F PW − F G |. The corresponding distributions are shown in Fig. S1, where in addition to the histograms of the module of the deviations it is also shown the histogram of the absolute force modules computed at the PW level |F PW |. Some statistical descriptors of all the distributions are listed in Table S1 to facilitate the ranking of the different basis sets in order of accuracy.
It is clearly seen that the force deviations for most of the Gaussian basis sets are very small compared to the absolute exact forces, being TZVP the only exception which performs clearly much worse than the rest. From the most to the less accurate set, they can be ranked as: On the other hand, the computational cost for each basis would rank them as TZV2PX-MOLOPT > TZV2P-MOLOPT > TZVP-MOLOPT ≫ QZV3P > TZV2P > DZVP-MOLOPT-SR > TZVP.
In view of these results, the chosen basis for the full study of the diffusion coefficients of H 2 O, H 3 O + and OH − was TZV2P which showed a very good balance between accuracy and computational cost.

II. EVOLUTION OF THE TEMPERATURE DURING THE SIMULATIONS
The dependence of the diffusion coefficients on the temperature of the system, in addition to the necessity of checking whether a proper sampling is being established in the 2genCP simulations, make it necessary to closely monitor the behavior of the temperature in all the cases. In Figs. S2, S3 and S4 I show the running average of the temperature ⟨T ⟩ and the evolution of its variance σ 2 T during the simulations for the neutral, acidic and basic water systems, respectively. The extent in which a proper canonical sampling is being achieved in a given NVT simulation can be checked by comparing the observed σ 2 T to its theoretical value in a NVT ensemble, i.e. σ 2 T (NVT) = 2 3N ⟨T ⟩ 2 [1]. This value for T = 300 K is represented in Figs. S2, S3 and S4 with a dashed black line.  Theor. NVT Figure S4. Evolution of the temperature in the simulations of the basic water system.