Explicit proton transfer in classical molecular dynamics simulations

We present Hydrogen Dynamics (HYDYN), a method that allows explicit proton transfer in classical force field molecular dynamics simulations at thermodynamic equilibrium. HYDYN reproduces the characteristic properties of the excess proton in water, from the special pair dance, to the continuous fluctuation between the limiting Eigen and Zundel complexes, and the water reorientation beyond the first solvation layer. Advantages of HYDYN with respect to existing methods are computational efficiency, microscopic reversibility, and easy parameterization for any force field. © 2014 Wiley Periodicals, Inc.


Introduction
Proton transfer is a crucial step in many chemical and biological processes, but difficult to observe experimentally at the atomic level. Instead, most knowledge about the mechanistic details of proton transfer are obtained from computer simulations and several methods have been developed to simulate proton transfer that vary greatly in their complexity, accuracy and efficiency.
At present, ab initio molecular dynamics (AIMD) simulations provide the most realistic description, but also have the highest computational costs. Therefore, the system sizes and timescales that can be reached are restricted to a few hundred atoms for a few 100 picoseconds. Despite these limitations, ab initio simulations have provided important atomistic insights into the transfer mechanism of hydrated protons, [1,2] which can be used for developing and parameterizing more approximate approaches. In an effort to gain efficiency, while preserving the detailed description of the excess proton, approaches based on semi empirical methods [3,4] and hybrid quantum/classical protocols have been proposed. [5] More towards a molecular description, but still capable of explicit proton transfer, are empirical valence bond (EVB) [6,7,8,9,10,11] and reactive force field [12,13,14,15,16] approaches. Typically, a classical molecular dynamics (MD) simulation is performed, during which proton transfer is evaluated in an effective potential energy landscape extracted from ab initio calculations of model transfer events. The accessible system size and timescale are larger than in ab initio calculations, but still an order of magnitude smaller than in classical force field MD.
Finally, the computationally most efficient methods are classical MD simulations that incorporate instantaneous proton transfer moves. The excess proton is covalently bound to a specific atom, but can be transferred to another atom by a MC move (Q-hop), [17,18] or after satisfying a number of triggers (Reactive Molecular Dynamics, RMD) [19].
Although these methods reproduce the kinetic picture of proton transfer quite well, they are not microscopically reversible and consequentially cannot produce correct thermodynamic ensembles.
Despite the availability of such wide range of methods for simulating proton transfer, we are not aware of any fully atomistic description that achieves (i) efficient sampling of a converged ensemble in a sufficiently large system, while (ii) maintaining thermodynamic equilibrium rigorously throughout the simulation.
In this work, we develop and test a new framework to describe proton transfer that meets these criteria. Our approach, which we call HYDYN (HYDrogen DYNamics), is computationally almost as efficient as the instantaneous transfer methods, but, in contrast, generates a correct thermodynamic ensemble. Efficient sampling is achieved by simulating within the classical approximation and using empirical force fields to describe all interactions. At regular intervals in a HYDYN simulation a proton acceptor is selected from among all possible acceptors around the current donor, using a Monte Carlo criterium that guarantees detailed balance. In between selection steps, the excess proton evolves on the free energy surface associated with proton transfer between the donor and the selected acceptor, using λ-dynamics. [20] A proton transfer step is considered successful if at the end of this period, the proton resides on the acceptor. After this period, the evolution is terminated and a new acceptor is selected from among the molecules nearest to the molecule that now carries the excess proton. In this way, the excess proton can visit every protonatable site in the system, mimicking the Grotthuss mechanism. Because the selection and evolution steps maintain detailed balance, HYDYN simulations yield a correct thermodynamic ensemble.
To test our method, we simulated an excess proton in small water clusters, bulk water and at a water/air interface, using both non-polarizable and polarizable force fields. Since the excess proton in water has been studied extensively in experiments and computations, there is a wealth of energetic, structural and dynamical data for validating the results of the HYDYN simulations.
The experimental and simulation data indicate a continuously changing structure [21,22,23,24,25,26,27,28] of the hydrated proton, fluctuating between an Eigen complex [22,23,24,29,30,26,31,32,27,33,28] and a Zundel complex (either as H 5 O + 2 or extended to H 13 O + 6 ) [30,26,34,35,24,31,27,28]. The limiting symmetric complexes are essentially never formed due to the asymmetry of the surrounding solvent [24,36,27]. Rather, intermediate asymmetric complexes dominate the structural ensemble, preferentially forming a distorted Eigen complex with one shorter hydronium-water hydrogen bond (called the special pair) [26,22,24,21,30,27,37]. The same complex can also be considered as an asymmetric Zundel complex where the central hydrogen is no longer centrosymmetric between the two oxygens [22,36,35]. The identity of the water involved in the special pair is continuously changing and performs what is called a "special pair dance" [22,23]. Proton transfer starts when the "special pair dance" ends, and a specific water persistently forms the shortest hydronium-water hydrogen bond. Then, the proton transfer takes place via a short lived Zundel complex to form a new Eigen complex with the former special pair water [25,22,38,23]. During proton transfer there is a collective reorganization of a larger number of water molecules and hydrogen bonds observed (up to 15 molecules) [39,32,40].
For the small water clusters, the interaction energies between the surrounding water and the Eigen or pseudo-Zundel complex calculated with HYDYN agree with ab initio quantum chemistry calculations. In addition, HYDYN reproduces the results of CPMD [21,22,23,24,25,26,27,28,40] with respect to the structural and dynamical properties of the excess proton in bulk water. In both CPMD [21,22,23,24,25,26,27,28,40] and HYDYN, the location of the excess proton fluctuates continuously between the limiting Eigen and Zundel states, as is the identity of the proton receiving site (special pair dance). Our bulk water simulations corroborate earlier observations that orientational changes of water molecules beyond the first solvation shell are important for the transfer process. Finally, the density distribution of the excess proton in a water slab obtained from HYDYN simulations yields surface affinities that are in good agreement with earlier studies [41,42,43].

Theory
To clarify the notation, we will first summarize the thermodynamic integration and λdynamics approaches. Subsequently, we will describe how we use λ-dynamics to model proton transfer between a donor and an acceptor in a classical MD simulation. Then, we will describe how the donor and acceptor are selected from among the possible donoracceptor pairs in the system. After presenting the basic concepts of HYDYN we will discuss the necessary parameterization and the technical details of our implementation in the Gromacs MD program. [44] Finally, we present a method for systematically analyzing the proton dynamics in a HYDYN trajectory.

λ-dynamics
Thermodynamic integration (TI) [45] is used to calculate the free energy difference (∆G) between a reactant state R and a product state P: Here, H TI is the Hamiltonian of the system, and λ is a coupling parameter that interpolates between the R (λ = 0) and P (λ) states, Figure 1: HYDYN concepts. (a) The two free-energy states of proton transfer between two water molecules. The proton can gradually change between fully interacting with the system (white) and non-interacting (grey) by slowly changing λ, thereby morphing between a hydronium and a water, respectively. Appropriate time evolution of λ can be achieved by applying λ-dynamics. (b) Proton transfer pairs considered during the selection procedure when the excess proton is located on the central water molecule.
To calculate ∆G via equation 1, λ is changed from 0 to 1 during the simulation, thus forcing the system from its reactant to its product state. The ensemble averages in equation 1 are then taken from the MD ensembles generated from the Hamiltonian H TI (λ) at different values of λ between 0 and 1. In the context of HYDYN, the reactant state defines the situation, in which the excess proton is bound to the donor, while the product state defines the situation, in which the excess proton is on the acceptor (Figure 1a).
Following the notation of Kong and Brooks, [20], we split the Hamiltonians of the reactant and product in λ-dependent (H R andH P ) and λ-independent (H env ) parts: In the λ-dynamics approach, [20] a Hamiltonian similar to equation 3 is used. In contrast to TI, λ is defined as an additional dynamic degree of freedom of the system with mass m λ , Figure 2: Projection of the angular coordinate θ onto λ, now defined as λ = 1 2 + 1 2 cos θ.
coordinate λ, and velocityλ. Accordingly, the Hamiltonian of the system is now expressed by with a force acting on λ, where V (λ) is the potential energy part of the Hamiltonian in equation 4: in equation 4, 1 2 m λλ 2 is the kinetic energy term associated with the λ "particle". The λdependent potential term U (λ) is introduced as a biasing potential that is parameterized such that the free energy landscape for the λ particle (equation 6) is consistent with both the experimental rate of proton transfer and the properties of the transfer intermediates (i.e., a Zundel-ion). The biasing potential is discussed in more detail below.

Constraining the λ-interval
Since only λ = 0 and λ = 1 represent physical states of the system, it is essential that the λ space is limited to the interval between the two physical states. In order to constrain the space sampled by λ, we switch to a new dynamic angular coordinate θ. By this modification, the actual dynamics takes place in θ space, and λ is redefined as the projection of θ on the abcissa ( Figure 2): The force acting on θ is with V the potential energy of the system, as defined in equation 6. As explained in Donnini et al., the conversion to θ introduces an entropic contribution to the effective free energy governing the λ-dynamics that stabilizes the end parts of the λ interval by a few kJ/mol. [46] We will show later that the λ-interval can be further reduced to [0, 1 2 ] without loss of generality. This additional restriction offers certain advantages in the implementation of HYDYN and in the analysis of a HYDYN trajectory.

λ-dynamics thermostat
The temperature of the λ particle is kept constant during a simulation by coupling the particle to an external heath bath. In HYDYN, we use the Andersen thermostat, [47] but other thermostat can easily be implemented as well.

Proton transfer with λ-dynamics
Now that we have summarized the key concepts of λ-dynamics, we can develop our protocol for incorporating the Grotthuss proton shuttling mechanism in classical MD simulations.
Since we want HYDYN to be compatible with all available force fields, most of which do not allow covalent bonds to be broken or formed, the proton transfer cannot be described in terms of bond breaking and formation. Instead, we introduce a proton transfer coordinate λ, which describes the proton transfer between the donor and acceptor within a selected donor-acceptor pair. How such pair is selected from all available acceptors will be described in detail later.
At λ = 0, the proton is bound to the donor, the reactant state, and is described by the potential energy functionṼ R . At λ = 1, the proton is bound to the acceptor, the product state, and is described by the potential energy functionṼ P . Figure 1a illustrates the concept. By changing λ from 0 to 1, the proton is effectively transferred, but, instead of being translocated, the proton appears (grows) at the acceptor, while it simultaneously disappears from the donor. In HYDYN the λ coordinate is assigned a mass m λ and allowed to change continuously between the reactant and product state. As in standard λdynamics, this change is driven by the potential energy surface, V (λ) (equation 6), which is linearly interpolated between the reactant state (donor protonated, acceptor deprotonated, Figure 1a, top) and the product state (donor deprotonated, acceptor protonated, Figure 1a, bottom). The equations of motion for the λ particle (equation 4) and the Cartesian coordinates of all atoms in the system are propagated simultaneously. [20] Thus, during the interval in which the pair is active, the proton can exchange continuously between the donor and acceptor.

Force field corrections
To correctly describe the free energy difference between the reactant and product states of a selected donor-acceptor pair, we need to include within the Hamiltonian the contributions to the free energy due to the breakage and formation of chemical bonds, which are not described by the force field. In HYDYN these contributions are described by an additional term V chem (λ) in equation 6, which will shift the donor-acceptor equilibrium by a certain free energy ∆G chem : To determine ∆G chem for the proton transfer reaction D-H + A → D + A-H, we consider the equilibrium between the reactant and product states. If we assume that the force field provides a sufficiently accurate description of the free energy contributions due the environment on the proton transfer reaction, the contribution ∆G chem that is due to the bond breaking and formation, is independent of the environment. With this assumption, ∆G chem can be determined by computing the difference between the free energy obtained from a force field simulation (∆G FF ) by applying equation 1, and the free energy (∆G ref. ) obtained from ab initio calculations or experiment, either in solvent or in vacuum.
Subsequently, the term V chem (λ) is obtained by fitting a polynomial function to the difference as a function of λ, as described by Donnini et al. [46] If proton transfer occurs between identical molecules, for example two water molecules, the contributions due to bond breakage and formation cancel and ∆G chem = 0. Since we will focus on the excess proton in water, there is no need to consider ∆G chem further in this work. In order to apply HYDYN to systems containing also a protein, a parameterization of V chem (λ) for each combination of amino acids as well as for all amino acid-water pairs is required. However, the parameterization procedure needs to be done only once for a force field, and will be undertaken in future work.

Biasing potential
In the standard λ-dynamics approach, a biasing potential (U (λ), equation 6) is used to will be discussed in detail in the Methods section.

2.2.3
Restricting the λ interval to [0, 1 2 ] In contrast to standard λ-dynamics, in which λ can sample the [0, 1] interval, we restrict the sampling of λ to the interval [0, 1 2 ]. We enforce this restriction by swapping the potentials V R andṼ P in equation 6 and simultaneously change λ to 1 − λ if λ passes 1 2 . In HYDYN the potentials are swapped by exchanging the force field parameters that define states R and P (see below). As follows from equation 6, the total potential energy is not affected by the swap. To conserve also kinetic energy and momentum of the λ-particle, we reverse the velocity of λ as well:λ = −λ. The advantage of reducing the λ space from [0, 1] to [0, 1 2 ] is that exchanging the force field parameters at λ = 1 2 , automatically transfers the mass of the excess proton from the donor to the acceptor, as we explain next. Furthermore, the reduction of λ space also simplifies bookkeeping of proton transfer events, because only λ = 0 corresponds to a physical protonation state, either at the donor or the acceptor.
In contrast, when using the complete [0, 1] interval, both λ = 0 and λ = 1 correspond to physical protonation states, which would complicate the analysis. Thus, while analyzing the location of the excess proton during a simulation, only frames with λ = 0 need to be considered.

Instantaneous mass transfer
Finally, we have to account for the transfer of one proton mass during proton transfer.
Following the same strategy as for the potentials and continuously change the masses of the atoms belonging to the donor/acceptor pair as a function of λ, is not possible, because the integration of the Cartesian equations of motion require a division by mass. Therefore, if the mass of the excess proton approaches zero while it disappears on either the donor or acceptor, numerical errors occur. To avoid such numerical problems, we transfer the proton masses instantaneously between the donor and acceptor when λ passes 1 2 . Thus, the mass transfer occurs simultaneously with the exchange of force field parameters.
To conserve kinetic energy and momentum after the mass transfer, we correct the velocities of the atoms, whose masses were changed: Here i enumerates the atoms, whose masses were changed, v n and v o are the velocities before and after mass transfer, respectively, C the correction factor and v c the velocity correction vector, which are defined as and The sum is over all atoms M for which the mass was changed, and m o and m n are the masses of those atoms before and after mass transfer. The derivation of these equations is given in the Appendix. Since imposing periodic boundary conditions violates conservation of angular momentum, we do not correct the change in angular momentum due to the mass transfer.
We note that the velocity redistribution following equation 11 is also time reversible, and therefore maintains microscopic reversibility. We have confirmed time reversibility in numerical simulations.

Transfer pair selection
The λ-dynamics based proton transfer protocol described in the previous section, allows a proton to move between two predefined sites, the donor and acceptor. We will call this pair of sites a transfer pair. Clearly, proton transfer within this transfer pair is not sufficient to describe a "free" excess proton that in principle can visit every protonatable site in the system. For the proton to visit every protonatable site, the transfer pair needs to change repeatedly during the simulation. Thus, to model the proton transfer process in a realistic manner, selection of a new transfer pair must (i) result in a physically-possible proton transfer path and (ii) be microscopically reversible.
We assume that proton transfer follows the Grotthuss mechanism, which has two important characteristics. First, the transferring proton can be any of the titratable protons on the donor molecule. For example, a hydronium ion has three equivalent protons, each of which can transfer to another protonatable acceptor, while a the histidine side chain has two titratable protons. Second, the excess proton can only transfer from one titratable site on the donor to an unoccupied titratable site on an acceptor via a hydrogen bond.
The selection of a new transfer pair must reflect these two features of the Grotthuss mechanism. We emphasize that during λ-dynamics only one transfer pair can be active per "free" excess proton, but that at selection steps all possible transfer pairs of the protonated molecule are considered. Therefore, the molecule that carries the "free" excess proton is determined first. Every titratable proton on this molecule is considered as part of a distinct transfer pair, following the first feature of the Grotthuss mechanism. To complete the selection, for each of these protons an acceptor is chosen by finding the acceptor site nearest to the proton, in accord with the second feature of the Grotthuss mechanism.
Finally, we select one pair from all available transfer pairs to become active until the next selection step.
With this selection procedure, every titratable proton on a protonated molecule is allowed to transfer over the course of the simulation. After successful proton transfer (see previous section) the identity of the donor and acceptor are interchanged (donor becomes acceptor and acceptor becomes donor), and the next transfer pair is selected from the new donor's protons. Sequential selection and transfer steps enable the "free" excess proton to visit all protonatable sites in the system in way that is consistent with the Grotthuss mechanism.
The selection of transfer pairs needs to conserve Boltzmann sampling, which can be achieved by applying an appropriate Monte Carlo criterium to accept or reject the new transfer pair. As was described above, the new transfer pair needs to include one of the available titratable protons on the protonated donor molecule. For each of the n titratable protons on a donor, the nearest proton-acceptor is selected to form a transfer pair. From these n pairs, one pair is selected at random with a normalized probability α i with E i the total energy of the system if transfer pair i is used to construct the λ-dynamics potential energy function (equation 6, Figure 1b). The sum in W tot runs over all transfer pairs available to the protonated donor, including the pair that has been active until the selection step. The pair selection is always enforced, but if the currently active transfer pair is selected, the system remains unaltered (equivalent to a rejection). Such a Monte Carlo move is equivalent to that used in all-exchange parallel tempering proposed by Calvo [48].
The probability to select a specific transfer pair i thus depends on the Thermodynamic equilibrium is conserved if the selection of the new transfer pair fulfills detailed balance. Therefore, we check that with N (o) the Boltzmann factor of the old pair that was active until the selection step By enforcing every selection, acc(o → n) = acc(n → o) = 1, and detailed balance is maintained. Including the previous transfer pair in α is essential to guarantee appropriate Boltzmann sampling, because this ensures that W tot is equal for the forward (o → n) and backward (n → o) move, and, hence that equation 17 is obeyed.
In contrast to instantaneous methods, such as Qhop or RMD, our transfer pair selection does not enforce an actual proton transfer, but rather a selection of which proton is eligible for transfer by λ-dynamics. This has two important advantages. First, the number and identity of the protons considered in the forward (o → n) and backward (n → o) move are always the same, maintaining detailed balance. In contrast, when actual proton transfer would be enforced an inequality could arise in the number of proton transfer steps considered in the forward and backward move. Such situation, which violates detailed balance, could occur, for instance, when protons migrate from bulk into a confinement, such as a pore. Second, selection of a new transfer pair involves a much smaller change of the system than a instantaneous proton transfer. Therefore, in HYDYN there is no need to minimize the energy (RMD) or to introduce a waiting interval after a move (Qhop), both of which also violate detailed balance.
During periods in which λ remains close to zero, the active transfer pair may switch frequently. We remark that such alternating pair selection in water resembles the special pair dance observed for the Eigen complex [26,22,24,21,30,27,37].

Technical details 2.4.1 Implementation into GROMACS
We have implemented HYDYN as described in the previous sections, in Gromacs 4.5.4.
We have used the existing free-energy routines and added the λ-dynamics functionality as well as the transfer pair selection procedure. To couple the λ-particle to a separate heat bath we implemented and applied the Andersen thermostat [47], which is suited for low-dimensional systems such as the λ-subsystem. [46]

Molecular topology
As explained previously, the potential energy surface on which λ evolves, is a linear interpolation between the potentials of the reactant and product states. In the reactant state the donor is protonated and the acceptor deprotonated, while in the product state the donor is deprotonated and the acceptor protonated. Therefore, we need two sets of parameters for each protonatable residue to describe the protonated and deprotonated state, respectively.
In Gromacs, both parameter set are included in one molecular topology and we kept that feature. For future extensions to aminoacids, which can bind protons at several distinct sites, such as Histidine, Aspartatic and Glutamic acid, more than two sets of parameters will be required.
For the excess proton in water, the two states represent a hydronium and a water, respectively. Therefore, each water/hydronium molecular topology consists of one oxygen with a charge q A (q B ) and van der Waals interactions σ A (σ B ) and  included.
By adding a suitable biasing potential, U (λ) in equation 6, we can bring the path closer to the physical path without affecting the thermodynamics. We achieve this by parameterizing the biasing potential U (λ), such that it (i) flattens the force field free energy profile (Figure 4), and (ii) introduces a barrier. To determine U (λ) we first computed the free energy of transferring a proton between two water molecules in water by means of the TI technique (equation 1, detail below) and fitted the following polynomial to the free energy profile with a, b and c fitting constants that flatten the force field's intrinsic free-energy profile ( Figure 4), and k a force constant that controls the Eigen to pseudo-Zundel ratio. We set k to 10 kJ mol −1 , which, if we define the Eigen state as λ < 0.25 and the pseudo-Zundel state as λ > 0.25, results in a free energy difference between Eigen and pseudo-Zundel of 2.7 kJ mol −1 , close to the ∼ 2.5 kJ mol −1 reported previously [22,21]. The barrier of 2.7 kJ mol −1 includes also the entropy effect associated with projecting λ on the angular coordinate θ (figure 2), which contributes 1.4 kJ mol −1 . [46] Note, that the total free-energy barrier along the λ-coordinate may not be directly related to the transition rate in the simulation, because other essential events, such as water reorientation, could be rate limiting instead.

Parameters that control proton transfer rate
After the biasing potential, U (λ), has been determined, there are two parameters, with which the user can control the rate of proton transfer in HYDYN simulations.
The first parameter is the mass of the λ-particle, m λ , which determines the dynamics Figure 5: Testing masses of the λ-particle and cut-offs of the transfer pair selection. a. The proton transfer frequency as a function of the mass and λ cut−of f . b. The evolution of λ in time with a mass of 0.001 a.u and a timestep of 2 fs displays a smooth time evolution not susceptible to integration errors. The results shown here were obtained with the SPCE water model [50] of the λ particle and thus the frequency of proton transfer. In this work, the mass was optimized such that it was sufficiently small to yield a reasonable transition frequency (Figure 5a), yet sufficiently high to use a 2 fs integration timestep (Figure 5b).
The second parameter controls the transfer rate via the pair selection. We introduce a function f (λ) that affects the probability to select transfer pair i (α i ) in equation 14, Since at the selection step f (λ) is the same for the forward and backward transition of the trial move, its choice does not affect detailed balance. However it does affect the probability to accept a trial move. Here we use a modified Heaviside function, or step function: where λ cut-off is the second parameter. The step function f (λ) ensures that a new transfer pair will only be selected when λ is smaller than the cut-off. The value of λ cut−of f is optimized based on two criteria. On the one hand, larger cut-off values lead to more successful proton transfer events (Figure 5a), because a new pair can be selected, even if the proton is not transferred completely onto the current donor (remember that because we restrict the λ interval to [0, 1 2 ], the situation in which the proton has remained at or returned to the donor, and the situation in which the proton has transferred to the acceptor are both characterized by λ = 0). On the other hand, to obtain a more realistic simulation, selection of a new pair transfer must only be attempted when the proton is (almost) completely located on one site forming an (almost) completely protonated molecule, i.e.
In our applications of an excess proton in water, we used a mass m λ of 0.001 a.u. and set the λ cut-off to 0.1. This choice was a compromise between transfer efficiency and proton localization ( Figure 5). If a larger timestep is desired, for example in combination with virtual sites, [49] the λ-mass has to be increased to prevent integration errors. The proton transfer rates can then be optimized by varying λ cut−of f ( Figure 5).

Analyzing a HYDYN trajectory
To analyze the simulations of an excess proton in water, we distinguish between the Eigen and Zundel complex. In our protocol, these complexes correspond to λ = 0 and λ = 0.5, respectively. To differentiate between Zundel and Eigen, we considered all frames with λ < 0.25 as the Eigen ensemble and all frames with λ > 0.25 as the Zundel ensemble.
In water, proton transfer events can be separated into complete transfer and rattle events [40]. We consider a transfer completed when three criteria are met: (1) λ must have passed λ = 0.5, which means that at least one parameter exchange and mass transfer should have occurred; (2) the identity of the molecule that is protonated has been changed; and (3) λ < λ cut−of f , i.e. the proton is (mainly) localized at one titratable site.
In a rattle event, the proton is rapidly moving back and forth between the same two protonatable sites. [40] A rattle event is present if two criteria are fulfilled; (1) λ must have passed λ = 0.5 more than once between two pair selection steps, thus at least two parameter exchanges must have taken place; and (2) λ > λ cut−of f from the first transition over λ = 0.5, i.e. the proton has not been fully localized at one of the two protonatable sites since it has left the donor.   [43], while the parameters for the TIP3P/SPCE hydronium were derived in this work, following the procedure outline in ref [43].

Methods
We used gromacs-4.5.4 [44] for the reference and a modified version of gromacs-4.5.4 for the HYDYN simulations. To model the excess proton in bulk water we placed one hydronium in a cubic box with 712 water molecules. To create a water slab, the z-component of this water box was increased to 20 nm. We used both fixed point charge and polarizable force fields, with the hydronium [43] dissolved in TIP3P [51], SPCE [50] and SWM4-NDP [52] water. The derivation of the non-polarizable and polarizable hydronium parameters in table i is described in a previous publication. [43] Bond distances were constrained using the SHAKE [53] and SETTLE [54] algorithm for hydronium and water molecules, respectively, and a timestep of 2 fs was used. The temperature of the system was coupled to a heat bath of 300 K, using the v-rescale thermostat by Bussi et al. [55] with τ t set to 0.5 ps. The λ-particle was coupled to an Andersen heat bath [47] of 300 K with a coupling constant of 0.2 ps. For the bulk water simulations the pressure was kept constant at 1 atm using the Berendsen barostat [56] with τ p set to 1.0 ps and a compressibility of 4.5 · 10 −5 bar −1 . For the water slab no pressure coupling was used. Van der Waals interactions were cut-off at 1.2 nm and the electrostatic interactions were treated using smooth PME [57,58] with a real space cut-off of 0.9 nm and a grid spacing of 0.12 nm. For the reference simulations, neighbor searching was performed every 5 steps with a 0.9 nm cut-off for the short-range neighbor list. For the HYDYN simulations, neighbor searching was performed every step.
The free energy profiles necessary for parameterizing the biasing potential ( Figure 4) were obtained by thermodynamic integration. We divided the interval in 21 evenly spaced λ points and at each λ point we sampled for 100 ns. The distance between the excess protons in the two involved molecule topologies was kept small by a distance restraints, which was zero until 1.5 nm, quadratic between 1.5 to 2.0 nm, and linear beyond 2.0 nm with a force constant of 1000 kJ mol −1 . The influence of this restraint on the free-energy profile was neglected, because the distance was normally within the 1.5 nm that constitutes the zero force region.
The water dissociation energy profiles were calculated following the procedure described in reference [43]. In short, starting from an ab initio optimized complex of a hydronium   [72]. We performed all ab initio calculations with the Gaussian03 package [73]. The force field energy profiles were calculated based on the B3LYP optimized structures.

Results & Discussion
To validate the HYDYN methodology, we considered three test systems. (1) A small cluster consisting of the Eigen or Zundel species plus the first water shell in vacuum.
The small system size allows us to perform energy calculations both with HYDYN and quantum chemistry approaches, so that we can validate HYDYN interaction energies. (2) The excess proton in bulk water. The availability of experimental and simulation data on this system enables us to validate the structural and dynamical properties of an HYDYN simulation. (3) The excess proton in a water slab. The excess proton is known to favor the interface over the bulk phase [41,42,43], which is an important thermodynamic property that HYDYN must accurately predict.

Energetics of the excess proton clusters in HYDYN
In the classical simulations that provided the basis for HYDYN, only the Eigen complex has been considered. Therefore, the hydronium parameters have been derived to best represent this complex. Since HYDYN also includes transfer intermediates such as the Zundel complex, accuracy of the interaction of this species with the waters in the cluster needs to be validated.
In HYDYN the real Zundel complex does not exist. Instead, proton transfer proceeds through a pseudo-Zundel complex that is the result of combining two hydronium topologies. The non-bonded interactions of the excess proton are not originating from one site centered between the two oxygens, but from two hydronium hydrogen sites that each interact with half their normal interaction strength. Therefore, the intramolecular properties To test the intermolecular interactions, i.e. the interactions of the pseudo-Zundel with the environment, we compare the energy profile of dissociation of the first-shell waters from the pseudo-Zundel core ion (H 5 O + 2 ) in vacuum and find good agreement between HYDYN and QM calculations, Figure 6. The deviation in energy is similar to that observed for the Eigen complex, Figure 6. Thus, our pseudo-Zundel species, which is basically a superposition of two hydroniums, interacts with the water environment in a manner very similar to a real Zundel species, in which the proton is shared between two water molecules.

Proton transfer in HYDYN
To verify whether HYDYN describes the properties of the excess proton in water, we performed HYDYN simulations of the hydrated proton in a small water box and carefully analyzed the trajectory. An example of a proton transfer event in our simulation is shown in Figure 7. As indicated by the degree of transparency, the excess proton disappears from the left and appears at the right water molecule, effectively transferring. In general, most proton transfers are not so smooth, displaying strong fluctuations in the λ-coordinate (see for example Figure 5b). This compares well to the picture of a continuously interchanging asymmetric complex [21,22,23,24,25,26,27,28,36].      We extract various dynamic properties for comparison to experiments and quantum chemical simulations, table ii. The Eigen/Zundel ratio and the transition rate are part of the HYDYN parameterization (k barrier and λ-mass, respectively), and therefore show good agreement with the reference values. Also the rattle frequency is in agreement with the reference value [40], despite the fact that in HYDYN only a subset of rattle events is possible. Finally, the hydronium diffusion constants are significantly increased with respect to the classical simulations, but do not reach the reference value. We speculate that, because in HYDYN the excess charge is distributed over at most two water molecules, in contrast to charge delocalization over all surrounding water molecules, the solvation shell is stronger and some flexibility in the proton transfer directionality could be lost and lead to a reduced diffusion constant.
Proton transfer has been explained in terms of breaking and forming specific water hydrogen bonds in the first and second solvation layer [40,39,24,22]. To investigate these hydrogen bonds, we analyze the radial distribution of hydrogens around specific oxygen(s), Before discussing the changes in the water hydrogen bond network upon proton transfer, we first discuss some key features of the individual rdfs. In the Eigen state, the peculiar first peak in the rdf of hydronium (magenta) at ∼ 1.7Å is due to the virtual site hydrogen at the accepting water. This peak is quite broad because at λ = 0 this site is not yet interacting. As λ increases, this site starts interacting, slowly introducing a hydrogen bond interaction that causes the peak to become narrower and shift to a smaller value.
Note that this partially existing hydrogen bond is not a true hydrogen bond, but a result of the way the excess proton is shared between two waters in HYDYN. No other hydrogen bonds are accepted by the hydronium in the Eigen state, as is evident from the plateau at 1Åin the cumulative rdf (N (r)). In addition, the hydronium first solvation shell waters (red) accept one strong hydrogen bond from the hydronium, responsible for the first peak at Changes in the water hydrogen bond network upon proton transfer are reflected in the difference in the (cumulative) rdf when moving along the λ-values. We observe a gradual transition of the rdfs as the proton transfer proceeds from Eigen to pseudo-Zundel (λ = 0 to λ = 0.5), where the rdfs of initially distinct oxygens (full color rdfs in Figure   9b) properly converge to a single rdf (black lines in Figure 9b) arising from equivalent oxygens. In addition, these graphs can also be interpreted in terms of a complete proton transfer. For example, in the upper graph, the hydronium oxygen gradually changes into an accepting water (magenta to black to red) or, vice versa, the accepting water changes into a hydronium (red to black to magenta). From the slight rise going from magenta to red in N (r) of the upper graph of Figure 9b we see that as the proton transfers to a neighboring water, the oxygen forms on average only slightly more hydrogen bonds.
However, the oxygens of the water molecules in the first solvation shell acquire on average almost a complete hydrogen bond, red to blue in middle N (r) graph in figure 9b.
Thus, upon successful proton transfer two of the waters in the second shell become first shell water (dark blue becomes red in Figure 9a), and, according to the rdfs depicted in the middle graph in Figure 9b, accept one hydrogen bond less. To accomodate this change, the hydrogen bond donating waters reorient. Simultaneously, two of the first shell waters become second shell water (light red becomes blue in Figure 9a) and accept an additional hydrogen bond, which again is accompanied by the reorientation of some waters. This observation is in very good agreement to other studies [40,39,24,22].

Equilibrium distributions of the excess proton
To demonstrate that HYDYN yields correct thermodynamic properties, we simulate an excess proton in a water slab. As shown in Figure 10, the probability density of the excess proton along the slab normal clearly shows an enhanced concentration near the surface.
The position of the maximum is comparable between HYDYN and classical simulations, but in HYDYN the affinity of the excess proton for the surface is higher. Differentiation intermediates explicitly into account. This is not surprising, considering the small free energy difference between the Eigen and Zundel states [22,21] and the resulting relatively large population of the transfer intermediates. A similar observation has previously been reported with respect to the proton hole [78].
The explanation for the higher surface affinity of the Zundel complex is likely the solvation structure in bulk. Both the Eigen and the Zundel complexes are poor hydrogen bond acceptors (see Figure 8 and 9), which strongly disrupt the water structure and enhance their surface affinity. The inability to accept a hydrogen bond affects two waters in the Zundel complex, as compared to one in the Eigen complex, and, hence, the driving force The enhanced diffusion of the excess proton also leads to faster convergence of the associated density distributions. The sampling time required to obtain the same standard deviation in the free energy difference between bulk and surface as in classical simulations is approximately half that of the classical simulation, as shown in Figure 11.

HYDYN in comparison to existing proton transfer methods
Since many approaches are available to model proton transfer events in MD simulations, we want to compare HYDYN with these methods. At the most detailed level there are ab initio approaches [79,1,2]. On a mechanistic level ab initio calculations are superior, but they come at a high computational cost, imposing harsh limitations on system size and simulation length. As a consequence, equilibrium sampling quickly suffers from convergence issues, making accurate assessment of thermodynamic properties difficult. HYDYN cannot compete on the mechanistic level, but the timescales and system sizes accessible with HYDYN are much larger.
Semi empirical methods, such as the approximate self-consistent charge density functional tight-binding (SCC-DFTB) method, [80,81] or methods based on the Modified Neglected of Differential Overlap, [82] like AM1 [83], PM3 [84,85] or OM n (n = 1, 2, 3), [86,87] often offer a good trade-off between computational cost and accuracy by simplifying the QM Hamiltonian. However, these methods do not provide accurate results for proton transfer reactions in water, primarily because they all underestimate the binding energy between water molecules. As a consequence, the first shell coordination numbers of the hydrated proton are too high, and the Zundel complex is overstabilized. [88,89,3]. Therefore, re-parameterization is essential, as was done by Wu et al. for the MNDO and OM n methods, [3] and Goyal et al. for the SCC-DFTB method. [4] However, because the computational costs of semi empirical methods remain much higher than of force field methods, the system size and timescales are more limited than in HYDYN simulations.
Reactive force fields [12,13,14,15,16] and EVB [6,7,8,9,10,11] approaches are computationally less demanding than ab initio simulations, but still describe proton transfer explicitly and include transfer intermediates. Nevertheless, a major drawback of these methods is that every unique excess proton topology requires parameterization [8] that is much more complicated than in HYDYN. In addition, MS-EVB is not microscopically reversible, due to molecules moving into and out of the EVB region in the course of a simulation [9]. To minimize the energy drift, the parameters have to be tuned carefully for every system [9]. Finally, the integration timesteps that are used in these methods are typically much shorter than in HYDYN. In comparison to HYDYN, these methods describe proton transfer events more realistically, but at a higher computational expense and a more complicated parameterization procedure for new systems.
The most efficient methods to simulate a mobile excess proton involve instantaneous transfer events. Examples are Qhop [17,18] and RMD [19]. Because Qhop parameters are available for all amino acids and several other molecules, using Qhop will save considerable time in setting up a simulation. At present this is a major advantage, but since the parameterization procedure in HYDYN is straightforward, as explained above, we expect that in due time parameters will become available for all relevant residues in all major force fields. A major disadvantage of these instantaneous transfer approaches is that they either suffer from poor proton transfer rates or do not obey detailed balance. Therefore, sampling the correct thermodynamic energy landscape is not guaranteed. Another drawback is that the characteristic proton transfer intermediates such as the Zundel complex are not taken into account. In general, HYDYN is a bit slower than the instantaneous transfer methods (∼ 1.5 fold if we assume the instantaneous transfer steps do not generate overhead), but, in contrast to qhop and RMD, includes transfer intermediates and rigorously samples at thermodynamic equilibrium.

Conclusion
We have described a proton transfer protocol for classical MD simulations (HYDYN) that conserves thermodynamic equilibrium, is computationally efficient and can reproduce the key aspects of proton transfer. The HYDYN code embedded in gromacs 4.5 is available as supporting information. Application of the HYDYN methodology is in progress to assess the diffusion properties of a proton on a lipid bilayer [90].

Acknowledgements
This work received financial support from the Volkswagen Foundation, Alexander von Humboldt foundation (MGW), and Academy of Finland (GG).

Appendix
Upon redistribution of the masses, the kinetic energy and momentum must be conserved.
Therefore, we adjusted the velocities of the atoms of which the mass is changed: where v is a velocity vector and the subscripts o and n indicate before and after the velocity redistribution that follows the mass transfer. The correction constant C and velocity correction vector v c are obtained by requiring that the total momentum p and the kinetic energy E kin are conserved.