Multi-State Reactive Molecular Dynamics Simulations of Proton Diffusion in Water Clusters and in Bulk.

The molecular mechanics with proton transfer (MMPT) force field is combined with multi state adiabatic reactive molecular dynamics (MS-ARMD) to describe proton transport in the condensed phase. Parametrization for small protonated water clusters based on electronic structure calculations at the MP2/6-311+G(2d,2p) level of theory and refinement by comparing with infrared spectra for protonated water tetramer yields a force field which faithfully describes minimum energy structures of small protonated water clusters. In protonated water clusters up to (H$_2$O)$_{100}$H$^+$ the proton hopping rate is around 100 hops/ns. Convergence of such rates is already found for $21 \leq n \leq 31$ and no further speedup in bulk water is found. This indicates that bulk-like behaviour requires solvation of a Zundel motif by $\sim 25$ water molecules which corresponds to the second solvation sphere. For smaller cluster sizes the number of available states, i.e. the number of proton acceptors, is too small and slows down proton transfer rates. The cluster simulations confirm that the excess proton is typically located on the surface. The free energy surface as a function of the weights of the two lowest states and a configurational parameter suggest that the ``special pair'' plays a central role in rapid proton transport. The barriers between this minimum energy structure and the Zundel and Eigen minima are sufficiently low ($\sim 1$ kcal/mol, consistent with recent experiments and commensurate with a hopping rate of $\sim 100$/ns or 1 every 10 ps) which lead to a highly dynamical environment. These findings are also consistent with recent experiments which find that Zundel-type hydration geometries are prevalent in bulk water.


Introduction
Proton transfer (PT) is ubiquitous in physical, chemical and biological processes, including PT in liquid water, 1-3 enzymatic catalysis 4 or protein-assisted proton transport in membranes. 5,6 One of the greatest challenges in understanding and characterizing proton translocation is to link kinetic or spectroscopic information with structural changes and the underlying energetics. In many cases, experiments provide thermally and/or structurally averaged data which is difficult to relate with atomistic or molecular aspects involving the motion of the atoms on the relevant time scales, although exceptions to this exist. 7 Spectroscopically, PT has been characterized for a range of systems. 3,[8][9][10][11][12][13][14] Proton transfer in bulk water and in clusters involves two processes: oscillatory proton transfer within a Zundel (Z) motif H 5 O + 2 and transfer from one Zundel motif to another one which is the so-called Grotthuss mechanism. 2 Gas phase spectroscopic work for small molecular frameworks has found an empirical relationship between the maximum of the infrared (IR) absorption band and the barrier height of PT reactions in combined investigations of computational simulations and experiments. 12,14 Such models can, however, not be applied to the situation in the bulk, such as proton transfer in water, because the Grotthuss mechanism involves transfer of charge through space.
In aqueous systems, the network of hydrogen bonds is a decisive factor in controlling charge transport. 15 Specifically the dynamics in the second solvation shell around the protonation site has been found to be essential for spatial proton transport. 16,17 Grotthuss transport itself is considered a picosecond process 1,18 although direct experimental measurement of this time scale is challenging. 19 The picosecond time scale makes ab initio molecular dynamics (AIMD) an attractive alternative. 1,20,21 Still, such approaches are usually limited to relatively short time scales (tens of picoseconds) 22 and small systems due to the appreciable computational cost involved. Nevertheless, depending on the level of theory, such ab initio MD simulations can provide important direct and molecular level insight into the energetics, dynamics and mechanisms underlying PT in solution.
To study proton translocation in (extended) condensed phase systems and on longer (nanosecond and longer) time scales, there has also been considerable interest in developing empirical methods capable of simulating Grotthuss proton transport. One such method is based on the empirical valence bond (EVB) method. 23,24 Following this, models such as two-state EVB (TS-EVB) 25 and multi state EVB (MS-EVB) 26 One of the conceptual challenges with EVB-based methods is the need to parametrize the off-diagonal elements which describe the coupling between the states. In the simplest case, the coupling is a scalar but quite often they depend on one or a few selected, geometrical coordinates ("reaction coordinates") which need to be intuitively chosen or motivated a posteriori. Another potentially challenging aspect in concrete applications is the definition of the valence bond states themselves. 32 Another method for studying proton transfer in molecular systems is Molecular Mechanics with Proton Transfer (MMPT). 33 In this approach, an accurate multi-dimensional potential energy surface (PES) for the proton transfer energetics between an acceptor and a donor atom is parametrized from ab initio reference data for generic motifs: symmetric single minimum (SSM; reference system is O 2 H + 5 ), symmetric double minimum (SDM; N 2 H + 7 ) and asymmetric double minimum (ADM; NH + 4 OH 2 ). The parametrization is based on Morse functions which allow hydrogen/proton transfer within a donor/acceptor motif.
Here, MMPT is combined with multi state adiabatic reactive MD (MS-ARMD) [34][35][36] to follow proton transport in the condensed phase. The method is referred to as MS-MMPT in the following. The approach taken in the present work is bottom-up in that first the energetics of the Zundel ion is treated within the MMPT framework. Next, the possibility for proton transfer involving a hydrated Zundel ion is provided by using MS-ARMD. Finally, the energy function is reparametrized by comparing with experimental data from small protonated water clusters. This model is then applied to larger clusters and proton transport in the bulk. is θ, which is the H *

Methods and Development
It is advantageous to introduce the unitless progression coordinate ρ = (r · cosθ − r min )/(R − 2r min ) 33,37 with r min = 0.8Å which maps ρ to the interval [1,0]. This choice for r min guarantees that the entire attractive region is included in the PES scan for all configurations (R, θ). The MMPT potential is where V 0 (ρ, R) is the isotropic part and k · d 2 (d = r · sinθ) is a harmonic approximation to the bending mode of H * A . The radial dependence of V 0 (ρ, R) is represented as a superposition of Morse potentials. Such parametrizations have been used in previous work 12,14,38 to investigate spectroscopic features of PT in small molecular systems.
In the standard MMPT force field, the total energy for one of the two possible resonance states I and II of the Zundel ion, see Figure 1, is the sum of all bonded terms (indicated by solid lines in Figure 1) plus the energy of the Because the PT motif is explicitly parametrized, no independent H 3 O + ion exists in the present framework. Within MMPT the energetics involving the transferring H * A is described by energies fitted to reference ab initio energies whereas the remaining (bonded and nonbonded terms) are those from an empirical force field, such as CHARMM. 39 This is achieved by using a geometry-dependent switching function such as All terms that are removed in one resonance structure and generated in another one are treated in this same manner.
In addition to switching on and off particular terms it is also necessary to modify certain internal (bonded) terms between the two resonance structures. Considering resonance struc- , l H 2 O eq ) are the force field parameters for an O-H bond in a hydronium ion and in water, respectively. In order to bring Eq. 4 into a conventional force field form, i.e. E bond = k eff · (l − l eff eq ) 2 , the expression is used with 38 These modifications are applied to bond and angle terms (see Table S1). Such mixed energy terms have been employed previously for spectroscopic studies. 12,38 To correctly describe the non-planar structure of Figure   1), written as where k d is the force constant and Ψ i,j is the dihedral angle.
The resonance structures (see Figure 1) also influence the non-bonded parameters of atoms where

Multi-state MMPT
Up to this point, transport of H * A (see Figure 1) is only possible within a predefined donoracceptor motif. In order to allow diffusion of charge in a condensed-phase context, additional states have to be introduced. This is done following the general concept of multi state reactive MD (MS-ARMD) which determines the energy and weights of relevant bonding patterns ("states") based on an empirical energy function 34,35,41 which is MMPT in the present case.
For a system of n water molecules and an excess proton, a "motif" (or state) in the present context is defined as an H 5 O + 2 ion and (n − 2) water molecules, see Figure 2.
Atoms O C to O F are potential H * acceptors from the first solvation shell (within the red background) and atoms O a to O d are those from the second shell to give rise for a second-shell motif such as For spatial transport of an excess charge, possible H 5 O + 2 motifs are determined for a given spatial arrangement x of the atoms where x contains the Cartesian coordinates of all atoms. This is done based on geometrical criteria as described below. The total potential energy of one such state j is Following MS-ARMD, the total potential energy for a given configuration x is then a linear combination of the energies of the m candidate states where In Eq. 10, w j (x) is the normalized weight 35,42 which has also been used for mixing double many-body expansion PESs. 43 However, other mixing functions are possible such as a tangent hyperbolic 34 The gradient of the total potential energy is then readily available from which is required for MD simulations.
From Equation 11 it is understood that the state with the lowest H 5 O + 2 energy among all candidate states is the state which contributes most to the total potential energy. Because the force field is reactive, it is possible for other low energy states, such as shown in Figure   2, to become the minimum energy state during dynamics simulations which eventually corresponds to a surface crossing between the current and a new H 5 O + 2 motif. This corresponds to Grotthuss proton transport.  Figure S1.
For computational efficiency it is advantageous to first determine V intra for all candidate states m as this is an inexpensive calculation. If the weight of state i is above a certain threshold (typically w i ≥ 10 −11 ; the effect of this choice on the results was tested) this state is regarded as a potential candidate state and the full potential energy is computed according to Eq. 9. Therefore, the total potential energy for a given nuclear configuration where m eff is the number of states retained according to the criterion for a minimal weight. were run in a similar fashion (without constraint). The solvent water molecules were modeled by a flexible SPC/fw water model 47 but with the MP2/6-311+G(2d,2p) equilibrium (l eq = 0.966Å, θ eq = 106.6 • ) bond length. Snapshots for analysis were collected every 5 fs.

MS-MMPT Simulations and Analysis
For determining the proton hopping rates, the time-course accumulation function, h(t), is recorded 29 Here, t is the simulation time and δt is the lag time δt = 5 fs between the frames. Initially, 3 Results

Parametrization of the MS-MMPT Force Field
The parametrization of the MS-MMPT force field was carried out in two steps. combines parametrization with respect to computed and experimentally measured reference data which has also been found to be beneficial for other parametrized PESs. 48 13 The parameters for the 3-dimensional reactive MMPT potential (Eq. 2) were determined by fitting to ab initio calculations at the MP2/6-311G++(2d,2p) level. 33,49 The fitting quality is shown in the inset of Figure 3a (R 2 = 0.999) and the parameters are given in Table S2.
Next, the total energy for the entire H 5 O + 2 ion was fitted to reference data at the same level of theory using two resonance structures as shown in Figure 1. For the water molecule the non-bonded parameters of the SPC model 50 were used. The progress of this fit is shown in Figure S2 and the resulting parameters are summarized in Table S3. The quality (R 2 = 0.88) of this fit from the room temperature simulations is shown in Figure 3a.    To determine the nonbonded parameters (partial charges and Lennard-Jones parameters) for [H 2 O] 3 H + 2500 reference structures were taken from an MD simulation at 300 K with calculated at the MP2/6-311++G(2d,2p) level of theory. Here, W 1 and W 2 are the two water molecules, respectively. All parameters from the Zundel-ion were retained and only the nonbonded parameters for the H 3 O + ion were adjusted. The quality of this fit is shown in Figure 3b.

Improved Parametrization Based on Infrared Spectroscopy
Spectra from finite-temperature simulations: In a next step the infrared (IR) spectra for protonated water clusters were determined using MS-MMPT-MP2 from finite-temperature MD simulations (50 K with ∆E = 10 kcal/mol). The IR spectrum was calculated from the Optimized cluster structures using MS-MMPT-IR are also reported in Table 1  Infrared spectra were also determined for larger clusters (H 2 O) n H + n = 11 and 21 ( Figure   5, right panels). For n = 11, the PT band is well described by the simulations (2600 ∼ 2900 cm −1 compared to 2400 ∼ 2900 cm −1 from experiment). The water bending mode peaks at 1600 cm −1 which agrees with experiment. Similar to the spectrum for n = 4 from MS-MMPT-IR simulations, the water stretching bands are between 3600 ∼ 3800 cm −1 . In the experimental spectrum, however, the bonded water-OH stretches are located between 3100 cm −1 to 3700 cm −1 . Specifically, the doublet near 3700 cm −1 was attributed to OH stretches of "acceptor donor" (AD) and "acceptor acceptor donor" (AAD) water molecules.
The sharp feature near 3600 cm −1 corresponds to "donor donor acceptor" (DDA) water molecules within the H-bonded network. 55 These assignments are based on power spectra of the respective internal coordinate which has already been used in earlier work on protonated water dimer. 49 For n = 21, similar spectral patterns for OH stretches and water HOH bends of water molecules were found compared to those of n = 11. The high-frequency part of the spectrum is realistically modelled whereas the region between 3200 and 3500 cm −1 does not show the broad distribution observed experimentally. 54 However, recent experiments 58 and simulations 59,60 suggest that this part of the spectrum is affected by more involved couplings between overtone bending (2ν b ) and OH-stretching vibrations.
There are a number of ways how the calculated IR spectra can be further improved or why the present simulations do not capture all features of the experimentally observed spectral signatures. First, all bonded interactions involving water molecules employ a harmonic en-ergy function for the OH-stretch and the OHO-bend. Secondly, it has been recently found that some spectral features may originate from a Fermi resonance between the HOH bending overtone and the OH-stretches as in water clusters and in liquid water. 59,60 This finding in charged water clusters is quite recent; 58 it has been previously believed that these features are primarily due to OH-stretch motions. 55 Capturing a Fermi resonance is usually not possible from classical MD simulations. 38 Rather, an effective excitonic Hamiltonian is used to describe such effects. 60 Table S5) MS-MMPT-MP2 force field calculations yield results close to those from the reference MP2/6-311+G(2d,2p) calculations (and therefore to CCSD(T)-F12/aVTZ) except for one of the bending vibrations which are described by harmonic potentials. 51 Overall, results from the MS-MMPT-IR force field (red symbols in Figure 6 Tables S2 and S3) and is therefore computationally considerably more effective.
In summary, the MS-MMPT-MP2 parametrization was refined by comparing with the IR spectrum for (H 2 O) 4 H + to improve the force constants of a few essential modes. Comparison with experiment for larger clusters is still reasonable for n = 11 but fails for n = 21 mainly because the water OH-stretch vibration is described with a conventional harmonic potential instead of an anharmonic (Morse) term which leads to a considerably larger density of states.

Simulations for [H 2 O] n H + Clusters for n = 6 to 100
With the improved MS-MMPT-IR parametrization MD simulations were carried out for small to medium-sized protonated water clusters with a particular focus on the proton transfer process itself. Each simulation was 10 ns in length and was run with both force fields at temperatures of ∼ 300 K using a range of ∆E values. All simulations were carried out in the N V E ensemble and the results reported below are averages over 10 independent simulations.   Table   S8.
For the smallest clusters the maximal transfer rates are smaller compared with larger clusters (see Figure 8) because the number of available proton accepting water molecules is small which also reduces the number of candidate states. It is found that for 21 ≤ n ≤ 31 the maximum rate, also representative for simulations in solution, is achieved. This suggests that "bulk-like behavior" is obtained for a finite cluster size of ∼25 water molecules hydrating the Zundel motif. This is consistent with earlier findings that the dynamics in the second solvation shell is essential for spatial proton transport. 16 For larger clusters the number m eff of retained candidate states saturates which leads to saturation of the hopping rate. From g H * Ow (r) in bulk simulations the number of water oxygens to form a cluster with n = 21 is found to include the second maximum and corresponds to the second solvation sphere. In other words, bulk behavior is found for solvent molecules including the second hydration shell.
The dependence of the hopping rate on the value of ∆E is related to the topography of the underlying potential and free energy surfaces. Therefore, G(q 1 , q 2 ) depending on two progression coordinates was constructed for different values of ∆E from unbiased MD simulations.
One progression coordinate (q 1 ) is the difference between the weights w 1 and w 2 of the lowest states as they distinguish between the Z-and E-states. For |w 1 − w 2 | = 1, the state is pre- Figure 8: Dependence of the PT hopping rates on ∆E from MD simulations using the MS-MMPT-IR force field for protonated water clusters with n = 10, 31 and 50. Data was averaged over 10 independent trajectories, each 1 ns in length. For these simulations and ∆E = 7 kcal/mol the maximum hopping rates are reached for all selected cluster sizes. For n = 10, the maximum rate is 69 hops/ns compared to 130 and 132 hops/ns for n = 31 and n = 50. Hence, the hopping rates depend on the cluster size but appear to converge with larger n.
dominantly that of a Zundel (Z) structure whereas for w 1 − w 2 = 0 (i.e. w 1 = w 2 ) there are at least two states that contribute equally, indicative of an Eigen (E) structure (see Figure   4). As a second progression coordinate q 2 the minimum δ of the displacement coordinate 20 and similar for δ AC and δ AD is considered. With these two coordinates, two-dimensional free energy surfaces (FESs) were generated from the probability For MS-MMPT-IR and ∆E = 5 kcal/mol only the Z-state is a minimum on the PES whereas for ∆E = 15 kcal/mol only the E-state appears, see Figure S3.   Table S8). The increase of the barrier for bulk water is reflected in the slowdown of the hopping rate to 60 to 70 ns −1 , depending on simulation conditions, see Table 3. Compared with the free energy surfaces for the two clusters that for the excess proton in bulk water exhibits an additional state, which is the "special pair", 68 discussed further below. Hence, proton diffusion and hops involve an additional state which may also contribute to the slowdown in going from the cluster to the bulk.

Proton Diffusion in Bulk Water
Next, the dynamics of the excess proton in bulk water was analyzed using MS-MMPT-IR.  The self-diffusion coefficient of the excess proton was determined from following the center of excess charge (CEC) 28,29 where w j is the weight of the j-th state and r j is the coordinate of weighted center of excess charge with q I/II k and r I/II k the charges and positions of atoms in the hydronium ion described in the two resonance structures. Figure 10 shows one trajectory following the CEC from a 250 ps simulation using MS-MMPT-IR which shows that extensive proton diffusion occurs. From r CEC the self-diffusion coefficient is obtained from its mean square displacement It is important to note here that in MS-MMPT no individual hydronium ion serves as a proton carrier. Rather, a Zundel motif as a whole carries the proton.
Hopping rates are again determined from time accumulation functions. Table 3 shows that the hopping rates in bulk simulations also depend on the value of ∆E, as was the case for the clusters (see Table 2). They range from 20 ns  Table S9) and are close to those for the n = 100 cluster simulations (i.e. 130 hops/ns). The diffusivity of the center of excess charge approaches a value of 0.30Å 2 ·ps −1 . Figure 10 reports the relationship between D CEC and the hopping rate. It can be seen that with increasing hop rate the diffusivity increases, too, but at a slower rate.
In a typical simulation, the average number of candidate states is 16.7 for first-shell motifs and increases to 110.0 if water molecules in the second shell are considered in the state selection (see Figure S5). Application of the weight criterion (w i ≤ 10 −11 ) to these states to retain them for computing the total potential energy of a configuration leaves an average number of 8.3 and 13.0 when including first-and second-shell solvent molecules, respectively. Including second-shell motifs is essential for conservation of total energy (see Figure S6) in particular when the identity of the minimum energy state changes. Increasing the time step to ∆t = 0.5 fs leads to an energy drift of 0.0025 (kcal/mol)/ns/molecule, compared with 0.0025 (kcal/mol)/ns/molecule from MS-EVB 3.2 simulations and 0.0028 (kcal/mol(/ns/molecule for aMS-EVB 3.2. 69 Table 3: Self-diffusion coefficients (Å 2 ·ps −1 ) from bulk simulations with two box sizes using the MS-MMPT-IR force field.
System size ∆E kcal/mol ∆t fs Similar to the cluster simulations, Figure 9 (bottom row) reports G(q 1 , q 2 ) for bulk simulations using MS-MMPT-IR. The FES exhibits three minima -the expected Z and E states and an additional third state for which the geometry (reflected in the value of δ) is that of E but the weights of the two lowest states are those of a Z state. This is the "special pair" (for a geometrical analysis of the special pair see Figure S8) known from previous MS-EVB found to be more stable than the Z-state, by between ∼ 0.5 and ∼ 1 kcal/mol, depending on the progression coordinate used. 69 The relative stabilization of and barrier between the "special pair" and the Z-state have also been determined from free energy simulations using the DFTB3-diag+gaus semiempirical method. With δ as the progression coordinate, it was found that the Z-state is destabilized by ∼ 0.4 kcal/mol and separated from the "special  73 notably the position of the most pronounced maxima except for the first maximum of g OwHw (r). The second maximum for g OwOw (r) is too faint and occurs at too large a distance r.
As a final experimental observable, the water-water and water-hydronium radial distribu-tion functions were determined (see Figure 11). First, the g OwOw (r) (water-water oxygen), g OwHw (r) (water oxygen-hydrogen) and g HwHw (r) (water-water hydrogen) radial distribution functions were obtained from 8 ns simulations in which the Zundel ion was removed.
The computed g(r) are compared with data from combined X-ray/neutron scattering experiments. 73 For g OwOw (r) the first density peak appears at r = 2.75Å which agrees well with experiment. The experimental data shows a second peak at 4.5Å, which is shifted to longer separations and less pronounced. For g OwHw (r) the first maximum occurs at r = 1.81 A from MD simulations, compared to 1.74Å from the experiment and the second peak at 3.28Å is well reproduced from both force fields. Similar agreement is found for g HwHw (r).
The RDFs were also calculated from simulations of the hydrated Zundel ion. From the 8 ns simulations the g O * Ow (r) (hydronium-water oxygen), g O * Hw (r) (hydronium oxygen-water hydrogen), g H * Hw (r) (hydronium-water hydrogen) and g H * Ow (r) (hydronium hydrogen-water oxygen) were determined and compared with experiment. 72 In all cases the position of the first peak agrees well with experiment. The small feature around 2Å in g O * Hw (r), not present in the experiments, has also been found in recent simulations using MS-EVB3. 2 69 and was related to the presence of weakly bound water molecules near the hydronium ion, which are also present here, see Figure S7. It should, however, be noted, that the MS-EVB 3.2 potential was modified to reproduce this feature The first minimum with MS-MMPT-IR and DFTB-diag+gaus is at r ≈ 4.0Å whereas that with MS-EVB is at r ≈ 3.5Å but this g(r) from experiment is unstructured in this region.
The difference between the calculated and the experimental g H * Ow (r) ∼ 2Å was found to also depend on the van der Waals radius of the O* (proton carrying oxygen) atom. Reducing this from 1.74Å to 1.65Å the maximum in g H * Ow (r) disappears and the shapes for g O * Ow (r) and g H * Hw (r) considerably changes, as demonstrated by Figure S9. Furthermore, the shape of g O * Ow (r) agrees better with experiment with smaller van der Waals radius of O*.
For comparing experiments and the present simulations it should be noted that in the experiments 72 the proton concentration is 50 times larger (1:10 for proton:water) than that used in the present simulations (∼ 1 : 500). Furthermore, experiments were carried out for concentrated HCl solutions, i.e. the counterions were present, and the experimental results (composite partial structure factors) need to be first processed by Monte Carlo simulations using an empirical potential structure refinement. 72 Hence, the differences between the computed and experimentally observed g(r) are not too surprising and it would, in fact be problematic if the g(r) from experiment and simulations did agree given the differences in system and data processing between the two approaches.

Discussion and Conclusions
The present work introduces a computational model to study proton diffusion in bulk water based on a bottom-up approach starting with extensive reference data from MP2 The present approach can also be compared with multi state EVB. 26,69 For the proton hopping rates the present simulations (100 to 130 ns −1 ) support results from MS-EVB3 simulations for which 108 hops/ns (see Figure 10) were found 29 whereas semiempirical DFTB, which underestimates the barrier for proton transfer, and CPMD simulations find between 400 and 500 hops/ns. 29 Given the entirely independent simulation technology and parametrization scheme the consistency of MS-MMPT-IR and MS-EVB3 lends support for a slower proton exchange rate (∼ 100 ns −1 ) than that found from ab initio (semiempirical or DFT) simulations. Based on the good agreement with available experimental data and existing results from MS-EVB3 we consider MS-MMPT-IR with ∆E = 7 kcal/mol to be the best model so far.
There are several possible improvements to the present parametrization. One concerns the parametrization of the solvent water molecules. Instead of harmonic OH-bonds an anhar-monic parametrization akin to a Kumagai, Kawamura and Yokokawa potential 76,77 could be considered for improved spectroscopic properties. Furthermore, recent work on the terahertz and Raman spectroscopy of pure water suggests that small amounts of charge transfer between H-bonded water molecules further improve the spectroscopy and dynamics of water. 57 Finally, due to the strong local electric fields around a Zundel ion, polarization effects are expected to be appreciable and including them should further improve the model. 7 The dynamics on the present MS-MMPT parametrizations was followed by using classical MD simulations. However, quantum effects have also been included for quite some time. 20,78,79 Using centroid MD simulations and the MS-EVB3.2 parametrization, one of the findings is that the self diffusion coefficient increases from ∼ 0.3Å 2 /ps to ∼ 0.5Å 2 /ps, which is an increase of about 70 %. This general finding is also consistent with recent ab initio ring polymer MD simulations including quantum effects, which report that including quantum effects generally increases the diffusivity of the charge defect. 80 However, this work also found that the actual diffusion coefficient depends sensitively on the density functional chosen and it was concluded that ".. it is difficult to gather enough statistics to estimate the diffusion coefficient from ab initio simulations using density functional theory, and when this is done the results turn out to be very sensitive to the choice of density functional." 80 Depending on the system considered, the diffusion coefficients ranged from 0.16 to 0.49Å 2 /ps. On the other hand, the influence of quantum effects on the radial distribution functions was found to be rather small. 81 Hence, depending on the property considered, including quantum effects may be required for semi-quantitative agreement with experiment.
With respect to simulation technology, EVB 23 requires diagonalization of a matrix with valence bond states as diagonal elements and off-diagonal elements that couple them at every time step. Conversely, MS-ARMD uses a linear combination of a number of candidate states.
The off-diagonal elements in EVB have led to some debate 82 as they need to be parametrized 37 in one way or another which often involves their dependence on geometrical coordinates. In MS-MMPT the energy mixing parameter ∆E is the only free parameter once the force field for the gas phase and hydrated Zundel ion are parametrized. From a conceptual viewpoint MS-MMPT mixes the energies of candidate PT motifs (see Eqs. 10 and 11). This differs from MS-EVB which is rather based on the concept of charge delocalization 69 which, in MS-MMPT, is a consequence rather than part of the model.
In terms of computational efficiency, MS-EVB requires one to solve an eigenvalue problem at every MD step and all states and coupling elements need to be computed. Contrary to that, in MS-MMPT a relatively large number of candidate states (∼ 100 in the present case) are considered which reduce to 10 to 20 on average for which the full force field energy needs to be determined. With regard to parametrization, MS-MMPT follows a bottom-up strategy whereas MS-EVB is rather more a top-down approach.
Given the considerable differences between MS-MMPT and MS-EVB3.2 it is gratifying to note that the two approaches arrive at similar conclusions, specifically concerning the hopping rates in bulk and the radial distribution functions although, specifically for g H * Hw (r), MS-EVB 3.2 compares better with experiment than MS-MMPT-IR. However, as already noted above, comparison with experiments needs to be done with circumspection as they use a higher proton concentration and counterion effects may lead to additional differences.
Further refinements (see Figure S9) of MS-MMPT-IR by optimizing the nonbonded parameters are possible but outside the scope of the present work. Also, the important role of the "special pair" is found with both approaches which lends credibility as to its relevance for a molecular-level description of the excess proton in water.
Recently, machine learning has emerged as a potentially interesting approach to represent high-dimensional (reactive) potential energy surfaces. 83,84 One of the hallmarks of such approaches is that they are very data-intensive. For a recent reactive MD study for the Diels-Alder reaction involving 3-dibromo-1,3-butadiene and maleic anhydride a NN based on ∼ 225000 structures was trained. 85 The evaluation time for this NN was a factor of ∼ 200 slower compared with a representation based on MS-ARMD which is the technique also used here. For a condensed phase system, such as the excess proton in bulk water where water molecules up to the second solvation shell (H 2 O) n H + with n ∼ 20) need to be included, the number of the electronic structure calculations will be considerably larger due to the flexibility and structural heterogeneity of the system. Hence, despite the potential benefits of such data-driven approaches, their use for reactions in solution, which involves even greater conformational variability than the Diels-Alder reaction, is still not routine.
Another advantage of MS-MMPT is the fact that it can be easily extended to chemically heterogeneous systems as the method is entirely based on site energies, i.e. the total energy of the system depends only on where the excess proton However, in the future, techniques such as transfer learning 86 may provide computationally advantageous routes to addressing this problem.
The analysis of the free energy surfaces underlines early findings which reported "...the com-plex behavior of the hydrated proton, bringing out features of both Eigen's and Zundel's views, but it also emphasizes that the defect cannot be fully understood using either of these...". 20 More recent spectroscopic work 7,66,67 reported that Zundel-type hydration geometries are prevalent in bulk water. This is consistent with the present findings using the MS-MMPT-IR parametrization which locate the "special pair", corresponding to a hydrated Zundel structure, as the minimum on the free energy surface, separated by 1.2 kcal/mol from the Zundel form itself. This is in good agreement with recent experiments which report a barrier for proton transfer of ∼ 1 kcal/mol. 87 Because the barrier between the two Zundel and the Eigen form is higher (1.4 kcal/mol) it is expected that the dynamics of the excess proton in bulk water is dominated by Zundel-like forms. This is also supported by recent experiments which found that the spectroscopic signatures for proton transfer are strongly IR active and weakly Raman active, indicative of an approximately centrosymmetric Zundel-like structure as the main spectroscopically responsive species. 67

Supporting Information
The supporting information contains detailed illustration to the MS-MMPT method and the MMPT potential for single proton transfer. Parameters of the MS-MMPT force field, normal mode analysis for small hydrated proton clusters and more details of MS-MMPT simulation for clusters and water bulks with one excess proton are also provided. All materials are available free of charge via the Internet at http://pubs.acs.org/.