Proton transfer in bulk water using the full adaptive QM/MM method: Integration of solute-and solvent-adaptive approaches
The quantum mechanical/molecular mechanical (QM/MM) method is a hybrid molecular simulation technique that increases the accessibility of local electronic structures of large systems. The technique combines the benefit of accuracy found in the QM method and that of cost efficiency found in the MM method. However, it is difficult to directly apply the QM/MM method to the dynamics of solution systems, particularly for proton transfer. As explained in the Grotthuss mechanism, proton transfer is a structural interconversion between hydronium ions and solvent water molecules.Hence, when the QM/MM method is applied, an adaptive treatment, namely on-the-fly revisions on molecular definitions, is required for both the solute and solvent. Although several solvent-adaptive methods have been proposed, a full adaptive framework, which 1 is an approach that also considers adaptation for solutes, remains untapped. In this paper, we propose a new numerical expression for the coordinates of the excess proton and its control algorithm. Furthermore, we confirm that this method can stably and accurately simulate proton transfer dynamics in bulk water.
Introduction<!>Theory and Method<!>Computational details<!>Results<!>Diffusion coeffcient<!>Discussions<!>Conclusions
<p>Proton transfer (PT) is an important phenomenon in biology, engineering, and solution chemistry. Although the first PT model was proposed 200 years ago by von Grotthuss, details of the mechanism of PT and solvation structure of the hydronium ion (whether it is an Eigen vs. Zundel cation) have not been fully understood because experimental observations cannot be straightforwardly interpreted, and such cations cannot be uniquely distinguished.</p><p>However, in the late 1990s and in the first decade of the 21st century, computational advancements have facilitated the development of plausible simulations and advanced knowledge of the mechanism. [1][2][3][4][5][6][7][8][9][10] PT is accompanied by the creation/annihilation of covalent bonds between oxygen and hydrogen, which is known as "structural diffusion," and which requires a quantum mechanical (QM) description that can explicitly count the electronic structures.</p><p>Therefore, ab initio molecular dynamics (AIMD) has played a major role in the progress of PT studies. In AIMD, the electronic structure of the entire system is evaluated to obtain the potential energy and forces acting on the respective atoms. 2,3,7,8,11 However, AIMD is severely limited in terms of the system size to be evaluated as well as the duration of the molecular dynamics (MD) simulation because of the significant computational cost. For instance, the relaxation of the solvation structure, including the second solvation shell, is assumed to be a rate-limiting process for PT. Therefore, the solvation shell should interact properly with the surrounding environment in the PT dynamics. Indeed, a recent study demonstrated that the radial distribution around hydroxide ions depends on the system size. 10 However, the extension of the system size results in a significant increase in computational cost. Obviously, large molecules such as proteins are also beyond the scope of AIMD, despite the high demand for them in fields such as chemical engineering and biochemistry. Car-Parrinello MD (CPMD),, 12 which is one of the most popular AIMDs, introduced a fictitious mass for electrons. However, an anomalously large value has been employed for the fictitious mass to ensure adiabaticity, which can distort hydrogen dynamics as deuterium. 9 It should be noted that divide-and-conquer (DC) treatment,, 13,14 in which the entire system is fragmented into a number of subsystems, reduces the CPU time (for instance, to O(n 1.2 )) when used in combination with the density functional tight-binding (DFTB) method [15][16][17] for homogeneous water systems, where n represents the number of water molecules. However, DC treatment still requires a large number of CPUs in accordance with the number of subsystems. In addition, DC treatment can cause discontinuities when particles diffuse across subsystem borders, although it is assumed to be mitigated to some extent by the introduction of buffer zones. Hence, there remain challenges for quantitative investigations and applications in the PT mechanism.</p><p>On the other hand, a hybrid simulation, namely the quantum mechanics/molecular mechanics (QM/MM) method, may be an appealing alternative because it can reduce the computational cost by partially applying QM calculations to the system. 18 Because the computational cost of the QM/MM method depends primarily on the size of the local QM region as opposed to that of the whole system, the QM/MM method has been widely employed in research on large molecular systems. However, the QM/MM method cannot be directly applied to PT dynamics simulations, requiring adaptive treatment for both the solute and solvent. In conventional QM/MM simulations for solutions, a solute of interest is placed at the center of the QM region so that it is surrounded by QM solvent molecules, and molecular definitions are fixed throughout the MD simulation. However, owing to free diffusion, the surrounding QM water molecules are replaced by other water molecules (MM).</p><p>Furthermore, the transferred proton itself is not consistent throughout the MD simulation because of structural diffusion, as explained by the Grotthuss mechanism. In such a case, the position of the hydronium ion can deviate from the QM center and diffuse across the of the excess proton and the surrounding solvent molecules should be adaptively updated during the MD simulation so that the hydronium ion is constantly located at the QM center surrounded by QM water molecules. Regarding the solvent diffusion problem, to date, several solvent-adaptive QM/MM methods have been proposed. [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] To understand solventadaptive QM/MM methods, it is useful to introduce the concept of QM/MM partitioning, which describes how the entire system is divided into QM and MM regions. In general, adaptive QM/MM methods are based on multi-partitioning approaches, where multiple QM/MM partitions are considered for every MD time step. While these partitions share the same QM solute molecule, they have different numbers and combinations of QM solvent molecules.</p><p>Potential energy and forces are independently evaluated for respective QM/MM partitions.</p><p>Finally, the resulting potentials or forces are linearly combined to evaluate the effective potential or forces that are used to update the coordinates for the MD simulation. In general, the adaptive QM/MM Hamiltonian can be represented as:</p><p>Here, r is the distance between particles, where the subscript ξ represents the QM center, and alphabetic subscripts represent particles. V (n) and σ (n) are respectively the potential energy and weight functions for the nth partitioning. The second term, which is called the "effective potential", is the weighted sum of the potential energies of all partitions. The third term is the "bookkeeping term", which was introduced to cancel out the artificial forces that arise from the derivatives of the weight function. 23,39 As a result, the effective force acting on the ith particle becomes</p><p>where</p><p>represents the force acting on the ith particle evaluated for the nth partitioning. Although adaptive QM/MM approaches enable the incorporation of quantum chemical 4 effects of solvation into the MD simulation, most approaches suffer from severe artifacts that are termed "temporal" and "spatial" discontinuities. Because we have discussed the discontinuities elsewhere, 40,41 we will briefly introduce them here. A temporal discontinuity is rephrased as a violation of the Hamiltonian conservation caused by discontinuities in the effective potential energy surface. Some solvent-adaptive QM/MM methods, such as the sorted adaptive partitioning (SAP) 21 and size-consistent multipartitioning (SCMP) 31 methods, are free from temporal discontinuity. However, spatial discontinuity is manifested as the monotonic drift of the bookkeeping term during the course of the MD simulation.</p><p>Although some ad hoc corrections have been proposed, 33,36,41 spatial discontinuities are inevitable for any QM/MM method because they arise from the unnatural manipulation of dividing a homogeneous solution into different layers. Therefore, for fairness, the static QM/MM method should also be subject to spatial discontinuities, and not only adaptive QM/MM. However, spatial discontinuities have not been extensively discussed because they are not considered a critical factor in static simulations. Note that compared to multi-size approaches such as SAP and difference-based adaptive solvation (DAS), 23 a spatial discontinuity can be suppressed by size-consistent treatment, in which the number of QM solvents is consistent among partitions. 41 Towards a solute-adaptive method, the first step is to numerically express the position of the excess proton, which is called an "excess proton indicator (EPI)." To date, several groups have addressed the development and application of EPI. [42][43][44][45][46] Previously, Chakrabarti et al. proposed the expression of the EPI ξ as</p><p>where r i is the coordinate of the ith oxygen atom. 44 The weight function for the ith oxygen atom can be written as</p><p>where r ij is the distance between the ith oxygen and the jth hydrogen atoms, and r o and d are parameters. The weight function W i defined by Eq. 4 indicates that the hydronium ion is more weighted in the EPI estimation in Eq. 3 than ordinary water molecules whose W i values are close to zero. Although this indicator has been used as the reaction coordinate in the post-processing of the MD trajectory for PT in membrane proteins, it cannot be used for an adaptive simulation of the bulk solution because of its discontinuity, scale-dependency, instability, and computational cost. The problems above mainly arise from the fact that the weights for ordinary water molecules distant from the hydronium ion in (3) are not necessarily zero. Note that the residual contribution is not negligible, and can cause critical errors when accumulated. Let us assume that a simulation cell is filled with water molecules, including a solvated hydronium ion. In the case of a small simulation cell, EPI in Eqs. 3 and 4 may be suitable, with EPI corresponding to the position of the hydronium ion throughout the simulation. However, as the system size increases, the indicator points to the geometric center of the simulation cell regardless of the position of the hydronium ion. This is because the accumulated residual contributions from ordinary water become comparable to that of the hydronium ion (scale-dependency). In addition, when the weighted water molecules move across the periodic boundary condition, the EPI becomes discontinuous, causing it to violate the Hamiltonian conservation. As a result, the temperature unnaturally increases, destabilizing MD simulations (discontinuity). Furthermore, even ordinary water distant from the hydronium ion can temporarily have a value as large as that of the hydronium ion when it forms strong hydrogen bonds. This results in an extraordinary displacement of the QM center, destabilizing the MD simulation (instability). It should also be noted that the computational cost to evaluate the weight of N water molecules is almost 2N 2 considering all possible oxygen-hydrogen pairs, and thus, its computational cost will become as large as the force evaluation, which is the bottleneck of MD simulations (computational cost). It is also notable that Pezeshki et al. formulated EPI as</p><p>where r D and r A j are the coordinates of the proton donor and the jth acceptor oxygen atoms, respectively. 47 Here, W is a normalization factor and W mj is a weight function of ρ mj , which is defined as</p><p>Here, H m represents the mth hydrogen atom bonded to the donor. Although the EPI is based on a continuous weight function W mj , the distinction between the donor and acceptor can cause discontinuities when the definitions are swapped, corresponding to PT. In addition, the EPI employs a threshold distance for acceptors to be considered, and hence becomes discontinuous when the acceptor list within the threshold is updated. Thus, when the EPI is used for MD control, e.g., the QM center position in the adaptive QM/MM method, the dynamics is not based on the continuous potential energy surface, resulting in a violation of Hamiltonian conservation. To avoid such problems, the proton indicator should consist of continuous functions of the coordinates of all particles, and it should also be able to detect local structural attributes around the hydronium ion. Although it is obvious that the molecular geometries of solvent water and the solvation structure in the vicinity of the hydronium ion differ from that of bulk water, it is not clear whether there exists any index represented by a continuous function (which can identify the position of the hydronium ion without noise and error throughout the MD simulation). Further, it is not clear how such an EPI can be used to control the dynamics of the QM center during the MD simulation. In most solvent adaptive QM/MM methods, the weight σ (n) in Eq.1 of each partitioning was evaluated based on the distance between the EPI and the respective QM and MM solvent water molecules. For the Hamiltonian conservation, ξ should be a variable of Hamiltonian Otherwise, there exists no Hamiltonian that derives the effective force in Eq. 2. Let U be the sum of the second and third terms in Eq. 1. For Hamiltonian conservation, the partial derivative of U should be exchangeable as</p><p>where a, b ∈ {x, y, z}. Therefore, the weight σ (n) in 1 should be a function of a coordinate for any particle that follows Hamiltonian dynamics. This condition is not a problem in the case of conventional solvent-adaptive QM/MM approaches because the QM center is fixed to a particle in the solute atom throughout the MD simulation, and σ is a function of the distance from atoms in other solvent molecules. In contrast, the EPI used in solute adaptive QM/MM is supposed to be a function of the coordinates of some particles, and thus there exists no conserved potential to satisfy Eq. 7.</p><p>To achieve an accurate and stable MD simulation, we propose a modified representation for the EPI. In addition, we propose a new protocol to control the indicator during the MD simulation, in which we introduce a virtual particle representing the QM center using constraint dynamics, which is called the RATTLE method. 48 Finally, we demonstrate the benchmark simulation for PT in bulk water, in which the Hamiltonian is well conserved, achieving a stable and durable MD simulation.</p><!><p>Excess proton indicator In the present study, we propose a new EPI ξ, which was used as the QM center in the SCMP simulation,</p><p>where r i and η i are the coordinates of an oxygen atom in hydronium ions or water molecules and their virtual sites, respectively. Thus, the EPI was defined as the weighted sum of internally dividing points between the oxygen coordinate r i and the virtual site η i over all solvent oxygen atoms, where W i is the weight of the ith oxygen atom. When the parameter c = 1, the EPI is a linear combination of virtual sites η i with weight W i . In contrast, when c = 0, the EPI is a linear combination of water oxygen coordinates r i . Even in this case, the EPI is indirectly subject to the water hydrogen coordinate through weight W i , as described below. Thus, when the value of c falls between 0 and 1, the EPI becomes a dividing point between the two positions (Figure 1). In contrast to the previous EPI in eq 6, the present EPI treats atoms without distinction between solute H 3 O + and solvent H 2 O.</p><p>The virtual site eta is introduced to stabilize the MD simulation in order to suppress the unordinary displacement, and details of the roles and effects are provided in the Discussion section. Here, we clarify the definition of the virtual site η. Let r iα be defined as the distance between the ith oxygen atom and the αth hydrogen atom, and throughout this paper, alphabetic and Greek subscripts represent oxygen and hydrogen atoms, respectively.</p><p>Then, its Gaussian-type score function φ (r) is defined as</p><p>where α and r 0 are the parameters (see the Discussion for details). Next, we assumed a normalized score φ iα as follows:</p><p>Then, the virtual site η i for the ith oxygen atom was defined using the hydrogen coordinates r λ as</p><p>Next, we define the EPI weight W in Eq. 8 and introduce a function ψ i for the ith oxygen by the summation of the spline function S 1 of the distance r iα between the ith oxygen and the αth hydrogen:</p><p>Here, S 1 satisfies the following boundary conditions:</p><p>where R small and R large are the bonding parameters that satisfy R small < R large . Here, S 1 (r iα )</p><p>is the bonding score between the ith oxygen and the αth hydrogen, which ranges from zero to one. Thus, ψ i can be regarded as a continuous expression of the number of covalent bonds of the ith oxygen atom (see Figure 1). Note that although we applied the spline curve here, the function S 1 is arbitrary unless it satisfies Eq. 13. Next, ψ i of the ith oxygen atom was used as an argument for another spline function S 2 (ψ i ), which monotonically increases between 2.0 and 3.0, and satisfies the following boundary conditions:</p><p>Based on S 2 , the weight W i for the ith oxygen atom is defined as</p><p>Because W i increases with the number of OH covalent bonds ψ, oxygen, and hydrogen atoms, the EPI mainly reflects the coordinates of atoms in the hydronium ion moiety. We found that S 2 could become zero for oxygen atoms except for several molecules located close to the hydronium ion by tuning parameters r 0 and α in Eq. 9 (see the Methods for detailed values).</p><p>This implies that when proper values are set to R small and R large in Eq. 13, the contribution to the EPI can be limited to only two or three oxygen atoms in addition to several hydrogen atoms close to the oxygen atom. This fact helps to exclude most solvent molecules from the EPI evaluation by saving computational cost, as discussed below.</p><p>Weight function in size-consistent multipartitioning (SCMP) method Here, we briefly review the SCMP method. The details of these functions are provided in the literature. 31,41,49 The weight function σ (n) in Eq. 1 is defined as below for the SCMP method.</p><p>Here, the EPI ξ was adopted as the atom in the QM center. O</p><p>where λ</p><p>Constraint for QM center Substituting the sum of the effective potential and the bookkeeping term into Eq. 7 and taking into account cancellation, the condition in Eq. 7 can be rewritten as</p><p>where a, b ∈ {x, y, z}. If ξ({r}) is a function of the other variables in the Hamiltonian, then</p><p>which does not satisfy condition 19. However, when ξ is an independent variable, Eq. 20 will be</p><p>Thus, σ will have symmetry with respect to the differentiation order. This indicates that the EPI ξ is represented by any particle that follows the Hamiltonian dynamics. Hence, we placed a dummy atom at the position corresponding to ξ. Note that the dummy atom does not directly interact with other particles but indirectly interacts through a constraint, as described below. We chose the mass of the dummy atom to be 1.0 × 10 −8 au, which is small enough to satisfy the constraint condition well, and does not affect other particles.</p><p>The velocity Verlet integrator was employed for MD in the present study, and thus, the RATTLE algorithm 48 was applied to control the constraint. In addition to Eq. 8, the constraint condition with respect to the velocity is given as</p><p>The MD algorithm for the ith particle is as follows:</p><p>where λ c and λ v are Lagrange multipliers determined from iterations to satisfy Eqs. 8 and 25, respectively The mass of the dummy atom was set to 1e-8 a.u., which is sufficiently smaller than other particles. Thus, this artificial treatment does not affect the dynamics.</p><p>We found that the use of a dummy atom smaller than the one set in the present study did not make any significant difference in the dynamics. However, a dummy atom with a larger mass may slow down the motion of the EPI and make the RATTLE iterations unstable, which can lead to the collapse of the MD simulation.</p><p>Adaptive Langevin thermostat We also carried out Langevin dynamics simulation to maintain a constant temperature, where the coupling strength with the thermostat adapts to the QM profile, 31,40 which is an index of the degree to which a solvent molecule behaves as a QM molecule. Using the velocity Verlet integrator for Langevin dynamics,, 50 the coordinate r and velocity v are propagated as</p><p>where m and γ are the mass and friction coefficients, respectively. F is a deterministic force derived from the potential function V , and R is the Gaussian random force, which is defined as:</p><p>where ζ is a random number that satisfies ζ = 0 and ζ 2 = 1. In the limit γ → 0, Eq. 26 reduces to the ordinary velocity Verlet algorithm for Hamiltonian dynamics. In the SCMP method, the QM profile ω i for the ith solvent molecule can be represented as:</p><p>where δ</p><p>behaves as a pure QM model, which can be rephrased as the ith solvent molecule, is defined as QM throughout all weighted partitioning. ω i = 0 when the ith solvent molecule behaves as a pure MM, which can be rephrased similarly. In the present study, we associated the QM profile and friction coefficient of the i-th solvent molecule as</p><p>where γ 0 is a parameter. Because ω i changes at every MD time step, γ i for the respective solvent molecule also changed. This enabled the MM solvent molecules to be fully linked to the thermostat, and thus, the coupling strength gradually attenuated when the solvent molecule approached the QM center.</p><!><p>The SCMP method was implemented in a local version of the GROMACS 5.0.7 package. [51][52][53] In all SCMP simulations, a total of 80 QM/MM partitions were considered, where the forces and energy calculations were carried out based on different partitions. The system consisted of one hydronium ion and 2047 water molecules in a periodic cubic box with a side length of 40.28 Å. Each partitioning has one QM solute hydronium ion and 32 QM solvent water molecules. In the SCMP, the transition parameters s QM , t QM , s MM , and t MM were set to 6.4, 8.4, 4.0, and 6.4 Å, respectively (see 31,40 for details). In the partitioning updating protocol, we allowed the updated partitioning to have a degree of order of 75% for efficiency, as detailed in our previous work. 40 We employed the SPC-Fw water model 54 for the MM water. For the QM part, we employed DFTB3 [15][16][17] implemented in GROMACS, as reported previously 55 using standard 3OB parameter sets. 17 The electrostatic interactions in the MM-MM and QM-MM models were calculated using the particle-mesh Ewald method. 56 Both electrostatic interactions and van der Waals interactions were damped to zero in the range of 8.5% to 9.0 Å. The electrostatic potentials on the QM atoms induced by the charges of the MM atoms were obtained using the smooth particle-mesh Ewald method with a switching function for electrostatic interactions (electrostatic embedding). After equilibration for several picoseconds (ps), all MD simulations were conducted for 100 ps with time steps of 0.25 and 0.50 fs for hydronium ion solution and bulk water simulations, respectively. For the control of the EPI, we chose r 0 = 1.3 Å and α = 0.129 Å for Eq. 9. In addition, we set R small = 1.20 Åand R large = 1.32 Å for (Eq.13) and ν small = 2.0. For temperature control, we employed the friction coefficients γ 0 = 100ps −1 and 10ps −1 for (Eq.29).</p><p>14</p><!><p>Energy and bookkeeping term For stable MD simulation and the production of an accurate ensemble, the MD simulation needs to conserve the Hamiltonian (total energy) throughout its entire course. As shown in Figure 2, we evaluated the Hamiltonian with different values of c in Eq. 8 under the microcanonical (NVE) condition in the present simulation. The Hamiltonian was well conserved in all simulations over 2.0 ps from the beginning in each simulation. Based on the fact that PT between the hydronium ion and water molecules was observed several times during this time period, we assume that the PT was simulated on a continuous energy surface for several ps. However, the simulations with c =0.2 and c=0.8 collapsed at 3.7 and 1.2 ps, respectively, because either the partitioning weight σ in Eq. 16 or EPI weights W in Eq. 8 became zero. However, although the simulations with c = 0.0 and 0.5 were sustainable for more than 4.0 ps, the energy conservation is violated after 2.0 ps in the respective simulations, which was caused by the extraordinary fast displacement observed of the QM center ξ. Note that these problems were never observed for the first few picoseconds in the simulations, but they commonly occurred in the next several picoseconds, which implies that some errors are implicitly accumulated during the course of the simulation time. To clarify this point, we evaluated the bookkeeping term in Eq. 1.</p><p>As we have reported in our previous study, 41 the drift of the bookkeeping term is related to spatial discontinuity. Notably, when c = 0.8, it drifted by 1700 kJ/mol for 1 ps under microcanonical conditions, which is approximately 200 times greater than the shift of the bookkeeping term for pure water or several monoatomic ion solutions. This discontinuity arises from the systematic difference in potential energies between QM/MM partitioning, whose QM regions are compact and fragmented. Therefore, we suppose that the large drift may reflect the dynamics of the hydronium ion because the diffusion of hydronium ions is much faster than that of water molecules and other ions, and the deformation of the QM regions occurs sooner. Taking into account the Hamiltonian conservation, the decrease in the bookkeeping term should correspond to the increase in internal energy, in line with the temperature increase observed under the microcanonical condition in Figure 2. To be precise, the released energy manifests as the acceleration of solvent molecules around the QM/MM border. Then, if the energy release is faster than its dissipation, the solvation structure and dynamics are significantly distorted, leading to the failure of the MD simulation. We found that the drift of the bookkeeping term became more distinct as the parameter c increased.</p><p>Indeed, MD simulations with larger values of c tend to become stuck earlier.</p><p>As we proposed,, 41 the drift of the bookkeeping term can be alleviated by introducing a correction potential U to cancel out the artificial diffusive force. However, the ad hoc approach is not promising for finding the optimal form, although it requires additional computational cost. As an alternative to both stabilizing the simulation and reproducing plausible PT dynamics, we employed the Langevin thermostat in an adaptive manner so that the coupling strength with solvent molecules gradually changes in accordance with the changes in molecular definition based on its distance from the QM center. To this end, we introduced the QM profile ω i as defined in (Eq.28), which indicates how much a solvent molecule behaves as a QM molecule. In the SCMP method, the QM profile averaged over the MD simulation smoothly shifted from 1 to 0 as the distance from the QM center increased, indicating that the molecular properties of solvent water gradually alternated between QM and MM. 31,40 Figure 4 shows the temperature over the course of the SCMP simulation with c = 0.0 employing an adaptive Langevin thermostat with the friction coefficient γ 0 = 100ps −1 . It can be observed that the adaptive thermostat controls the temperature well, maintaining it at the reference temperature of 300 K, enabling the MD simulation to be durable over several hundreds of picoseconds. It should also be noted that the spatial discontinuity distorted the dynamics of the solvent molecules located near the QM/MM border rather than the QM center. 41 The thermostat can also alter the dynamics, but under the present adaptive usage, the solute and solvent molecules in the vicinity of the QM center are free from the Langevin thermostat. Therefore, it is plausible to evaluate the dynamical properties from the trajectories obtained in the present simulation, as discussed below. As will be discussed later, the friction coefficient γ 0 does not affect the results of the hydronium ion simulations. Thus, most analyses in the present study are based on simulations with c = 0.0 and γ 0 = 10 ps −1 , unless otherwise stated.</p><p>Computational cost Figure 3 shows the computational time for one MD step of the full adaptive QM/MM method compared with the solvent adaptive QM/MM method, both of which are based on the SCMP framework. Here, the solute molecule is H 3 O + in the full adaptive QM/MM and H 2 O in the solvent adaptive QM/MM. Note that although the system size increases twice, the entire computational time remains at a comparable scale. This</p><p>shows the advantage of the QM/MM-based approach over QM (ab initio) MD simulation, which faces a drastic increase in computational time with system size. Figure 3 indicates that force and energy calculations, which include QM(DFTB3), MM, and QM/MM interactions amounts to ca. 80%. Because the SCMP method is based on all partitioning calculations, the computation cost for this part linearly increases according to the number of partitionings.</p><p>Note, however, that these parts of the calculation on respective partitionings are independent of each other. Hence, it is suitable for parallel computation. In addition, because all partitionings in the SCMP method have a QM region of the same size, it shows higher scalability in parallelization than the other solvent-adaptive QM/MM. As a result, the MD performance of the proposed method is comparable to that of the conventional QM/MM MD simulation. The size dependency of this part mainly arises from the MM part calculation, which is somewhat moderate compared to that of the QM calculation. The second bottleneck is the evaluation of the SCMP weight function σ in Eq. 16, which is shared by both solvent-adaptive and full-adaptive QM/MM. Note that the computational cost for this process can be independent of the system size because the progress function λ in Eq. 17</p><p>becomes constant when the solvent distance from the QM center is larger than the threshold, and thus it is sufficient to count a limited number of solvent molecules to evaluate the SCMP weight σ. The residual cost contains the EPI evaluation and RATTLE calculation, which are specific to the full adaptive QM/MM method. Notably, these specific parts are so small that the computation costs did not vary significantly between the solvent and full adaptive QM/MM methods.</p><p>Solvation structure Hereafter, O* denotes the oxygen atom nearest to the QM center.</p><p>While O* represents the oxygen atom in solute H 2 O in the ordinary solvent adaptive method, in the case of hydronium ion simulation, O* is assumed to be a part of the hydronium ion. around the hydronium ion shifted from 2.8 Å to 2.6 Å, while the first peak in the experiments was observed at 2.5 Å. 57 In contrast, DFTB2 simulations for hydronium ions showed bimodal peaks at approximately 2.4 Å and 2.8 Å, indicating that a Zundel-type structure was the dominant component. 58 On the other hand, the present O*-O RDF showed a distinct single peak, which is consistent with the DFTB3 simulations. 59 According to the previous AIMD simulation,, 2 such a result implies that an Eigen-type structure is the major component.</p><p>The height of the first peak of oxygen atoms in bulk water was 4.1, which was significantly larger than the experimental value of 2.5. 60 The first peak of the O*-O RDF is broader than the empirical potential structural refinement (EPSR) of the experimental data. 57 As a result, the coordination number was estimated to be as large as 5.4 by integrating the first solvation shell. As DFTB3, which was employed for the QM region, causes oversolvation and high density for bulk water, such properties appear to be carried over to the hydronium ion simulation. With respect to the O*-H RDF, the first peak reflects the hydrogen atoms that are covalently bonded to the hydronium oxygen. While the bulk water simulation showed a second peak of O*-H RDF at 1.8 Å, it was not observed in the hydronium ion simulations, indicating that the hydronium ion did not accept hydrogen bonds. Note that the previous AIMD simulation showed that the O*-H RDF of the Eigen-type cation had an additional peak at approximately 1.6 Å. 7 Although an Eigen-type cation was more probable in the present simulation, as shown in the next section, we could not find any additional peak in that distance range.</p><p>Energetics and dynamics Figure 6 shows the potential of mean force for PT compared with the potential energy in the gaseous phase, which is projected on two reaction coordinates, the distance between two oxygen atoms R O*O' and hydrogen displacement δ.</p><p>The hydrogen displacement δ is defined as the difference between the two distances as and friction coefficient γ 0 did not affect the potential mean force (see Supporting Information). While the Zundel cation was more stable than the Eigen cation in vacuo, the balance was inverted in the aqueous phase, which is in agreement with previous studies. 8,58 Figure 6 gives a clear picture of the PT mechanism in the Eigen to Zundel to Eigen sequence, where the energy barrier for the PT disappears as R O*O' becomes smaller than 2.43 Å. The present SCMP simulation estimated the energy barrier for PT to be approximately 2.7 kJ/mol, which was in good agreement with the DFTB3/3OB simulation. 59 Although DFTB3 leads to the overbinding of OH covalent bonds and the underestimation of hydrogen bonds as previously reported,, 17,41 the estimated barrier was in agreement with previous studies of AIMD with DFT using BLYP 7 and HCTH functionals,; 61 however, it was significantly smaller than those of MS-EVB2 and 3, which were estimated to be 8.4 kJ/mol. 62 It should be noted that the nuclear quantum effect lowered the PT barrier. 7,9,63,64 Moreover, the energy barrier of approximately 1 K B T at room temperature is supposed to disappear by the incorporation of the nuclear quantum effect, causing the topological defect to delocalize and transfer at a rate faster than 100 fs. 7,9 In contrast, the present profile showed a distinct free energy minimum corresponding to the Eigen structure observed during the resting state of PT. Figure 7 shows the time evolution of PT projected on two reaction coordinates. Here, the PT event was centered around a moment with δ = 0, and the trajectories were averaged over more than 3000 PT events. Note that the PT appeared as an event within 100 fs. It should also be noted that concerted oscillations were observed along the two reaction coordinates with a period of approximately 15 fs, where R O*O' decreases and |δ| increases are clearly coupled.</p><!><p>In general, a thermostat can affect dynamical properties. To reproduce plausible dynamics, we employed the adaptive Langevin thermostat, where the coupling strength to the thermostat of the respective solvent molecule is adaptively updated according to the distance from the QM center. To verify the influence of the thermostat, we benchmarked the diffusion coefficient of the QM water molecules in the MM water system in the conventional solvent adaptive QM/MM method by employing two approaches: linear fitting of the mean-square displacement (MSD) and the integration of the velocity autocorrelation function (VAC) of the EPI. As shown in Table 1, when γ 0 = 10ps −1 , the diffusion coefficient of H 2 O shows good agreement with the NVE and QM-MD simulations with the DFTB3 model, while it is underestimated when γ 0 = 100ps −1 . This implies that the dynamics with γ 0 < 10ps −1 are almost free from artifacts from the Langevin thermostat.</p><p>Next, we calculated the diffusion coefficient of hydronium ions with full adaptive QM/MM using MSD and VAC, where the EPI of the present study was used to estimate the position and velocity of the hydronium ion. For evaluation, we employed six and fifteen independent 100-ps-trajectories with γ 0 = 100 and 10 ps −1 , respectively. We emphasize that the availability of VAC analysis indicates the advantage of the present full adaptive QM/MM simulation, which makes the velocity of the hydronium ion directly accessible. The resulting values obtained by the two analyses showed agreement and smaller statistical error, which implies that sufficient sampling was achieved by the stable MD framework.</p><p>The obtained diffusion coefficients varied according to the friction coefficient, where 0.58 Å2 /ps with γ 0 = 100ps −1 , and 1.0 Å2 /ps with γ 0 = 10ps −1 . Considering H 2 O, the result with γ 0 = 10ps −1 is more plausible and in good agreement with the experimental value of 0.94 Å2 /ps, although we suppose that there appears to be error cancellation. This agreement is in stark contrast with previous studies, where there can be found a systematic underestimation as 0.40 and 0.36 Å2/ps by MS-EVB2,  61 These underestimations may be attributed to the limited size effect in the DFTB3-diag simulation. 62,65 On the other hand, the limited size effect is less likely in the present simulation because of the sufficiently large system consisting of 2048 H 2 O molecules employed by taking advantage of the QM/MM method. It is also insightful to compare the previous DFTB studies, although they employed a modified parameter 3OB(diag) for the OH repulsive potential, while 3OB</p><p>was used in the present study. In this study, the obtained values are larger than 0.66 ± 0.20 Å2 /ps) in the previous paper based on QM-MD simulation using DFTB3-diag,, 58 while it is in good agreement with DC-DFTB3 (diag). The disagreement between the DFTB3-diag and DC-DFTB3-diag was also expected to result from the limited size effect in the DFTB3diag simulation. However, this discrepancy cannot be accounted for by the size-limiting effect alone because the DFTB3-diag for the 128 water system resulted in a H 2 O diffusion coefficient of 0.38 Å2 /ps, which was larger than the diffusion coefficient of 0.19 Å2 /ps for DC-DFTB3-diag with 513 water molecules (Table 1). The deviation may also have resulted from either/both the protocol of the diffusion coefficient calculation and/or DC treatment, such as the discontinuity caused by particles crossing the subsystems. In a previous study, the hydronium diffusion coefficient was indirectly evaluated by the summation of the vehicular and Grotthuss diffusion coefficients, where the former corresponds to the water diffusion coefficient, and the latter is estimated using the PT pitch and rate, and the two dynamics are assumed to be independent. On the other hand, the present approach with the use of EPI allowed the direct evaluation of diffusion without the assumption of independence. In addition, DFTB3/3OB and DC-DFTB3-diag showed different diffusion coefficients of H 2 O, while they had a comparable PT barrier in free energy (Table 1) Hence, the agreement between the present result and DC-DFTB3-diag implies that the correlation of solvent-PT dynamics and/or the other factor of the rate-limiting process of proton diffusion in this model.</p><p>Previous simulations have proposed that the rate-limiting step is the change in the coordination number of water molecules in the solvation shell. 2,9,66 In addition, water reorientation related to molecular rotation occurred together with the formation/deformation of hydrogen bonds. 67 Table 2 shows the resulting values obtained by the explicit integration of the second rank auto-correlation function of the OH bond orientation of H 2 O using the solvent adaptive QM/MM (SCMP) method. Although the result appears to be free from thermostat artifacts, it showed that DFTB3/3OB significantly underestimated the relaxation time when compared to the experimental value. This result is in agreement with the fact that the underestimation of the hydrogen bond energy by DFTB3 results in a faster relaxation of the bond orientation of water compared with that in the experiments. 59 Consistently, the diffusion coefficient of DFTB3 water was significantly overestimated. 59,68 Therefore, we assume that the overestimation of the free energy barrier and the underestimation of the orientational relaxation canceled out, which is in agreement with experiments regarding the diffusion coefficient of the hydronium ion. Because the nuclear quantum effect is assumed to lower the free energy barrier, the balance of error will be altered when it is incorporated.</p><p>However, it is difficult to precisely predict that it can indirectly have competing effects on water diffusion with respect to surrounding water molecules. 69</p><!><p>To achieve a solute-adaptive QM/MM, it is essential to numerically express the hydronium position. The indicator should be capable of distinguishing the solvation structure specific to hydronium ions from ordinary bulk water. Intuitively, the number of covalent bonds appears to be such a structural measure, but it cannot be directly used as an EPI because it is discretized. On the other hand, ψ found in Eq. 12 can be regarded as a continuous expression of the number of covalent bonds of an oxygen atom. where the subscript represents its positional order from the QM center. With respect to the nearest neighboring oxygen, which is assumed to be a part of the hydronium ion, ψ 1 almost took the value of 3.0, although it was observed to be below 3.0, in agreement with Figure 8(a). The second nearest oxygen to the QM center, which forms the nearest neighboring water, showed a bimodal distribution with respect to ψ 2 . The first majority was at ψ 2 = 2.0, which occupied 88%, while the second majority appeared at 3.0, which amounted to 4%.</p><p>The ψ 2 value of 2.0 corresponds to the resting state of PT, where the hydronium ion takes the Eigen form. In contrast, ψ 3 equaled 2.0, most of the time, and it rarely became larger.</p><p>It is obvious that oxygen atoms farther away from the QM center than the third nearest oxygen atom have a ψ i value of only 2.0. Because of Eq. 14, the EPI weight W i in Eq. 15 of oxygen with ψ = 0 becomes zero, which indicates that the EPI was insensitive to distant water molecules from the hydronium ion. This is the most important attribute of an EPI to avoid scale-dependency and instability, as described above. Furthermore, this also helps to drastically reduce the computational cost of evaluating the proton indicators because the number of water molecules that need to be considered is very small. While the EPI in previous studies was a function of oxygen coordinates as defined in Eq.</p><p>3, the EPI in the present study was expressed as a linear combination of dividing points between water oxygen and virtual site coordinates η, as defined in Eq. 8. Note that when r = r 0 in Eq. 9, the Gaussian-type function φ has the maximum value. Thus, r 0 was tuned to be larger than the ordinary OH covalent bond length. As a result, as shown in Eq. 11, η i mainly reflects the hydrogen position that has the longest covalent bond with the ith oxygen atom. If the oxygen coordinates are used in Eq. 8, that is, c = 0, the EPI during the resting state of PT (which constitutes most of the MD simulation) would correspond to the oxygen in the hydronium ion. Hence, the EPI is displaced between two oxygen atoms at an anomalously high speed when the proton transfers. However, the EPI with c > 0 fluctuated even during the resting state of PT, reflecting the stretching of OH covalent bonds. On the other hand, when the proton is transferred, it is displaced at a moderate speed, which likely reflects the speed of the hydrogen transfer because the EPI was mainly dominated by a longest-covalently bonded hydrogen atom, that is, the transferred proton. Figure 9 shows the velocity distribution of the EPI sampled in course of MD simulation with c = 0.0 and 0.2.</p><p>Indeed, when c = 0.2, the probability of a large velocity is smaller than that with c = 0.0, although the probability of the middle-range velocity increased. Because anomalously rapid displacement of the EPI can destabilize the simulation causing discontinuity, the non-zero value of c appears to be more suitable. However, this contradicts Figure 2(a), where c = 0.0 showed the best stability and sustainability. To account for this, Figure 2(b-c) is insightful, where c = 0.0 resulted in the smallest drifts of the bookkeeping term. Therefore, small but constant fluctuations of the EPI coordinates can frequently cause the deformation of the QM region and the corresponding partitioning update, which accelerates the drift of the bookkeeping term and distorts the dynamics. The distorted dynamics can also induce anomalously rapid displacements of the EPI, which is more critical than that observed when c = 0.</p><!><p>In this study, we proposed a numerical expression for the excess proton indicator for hydronium ion simulations in bulk water, and implemented it in the SCMP code. Based on the results, we achieved a QM/MM simulation that was both solute-and solvent-adaptive, which we refer to as a "full adaptive QM/MM method", and we successfully demonstrated its stability and efficiency. Based on this new framework, we were able to confirm the total energy conservation of the PT simulation, that is, temporal continuity, with Hamiltonian dynamics under microcanonical conditions. However, the bookkeeping term and the tem-perature constantly varied owing to spatial discontinuities, for which we next conducted a Langevin dynamics simulation by employing the thermostat in an adaptive manner to retain plausible dynamics of the PT. We emphasize that the present computational approach is advantageous for evaluating the dynamical properties of PT in bulk water for the following three reasons. First, the position of the excess proton is defined with numerical stability and obtained on-the-fly through MD simulations. This provides direct access to various physical properties, such as the diffusion coefficient through the MSD or velocity autocorrelation function, which have been indirectly evaluated in previous studies. Moreover, the EPI may be used as a reaction coordinate for enhanced sampling, such as umbrella sampling, by imposing an artificial force on a dummy atom. Notably, the present EPI can also be used for post-MD analysis, as long as the coordinate information is retained. This may be useful to reassess previous AIMD simulations by removing artificial noise. Second, the computational cost required for the present approach based on the QM/MM method is moderate compared to AIMD, making it accessible to longer dynamic trajectories. As a result, statistical errors can be significantly suppressed compared to conventional studies using the AIMD method. The linear shape of the obtained MSD and the agreement between diffusion coefficients obtained by MSD and velocity autocorrelation function provide strong evidence that statistical errors are greatly suppressed. We assume that this feature will further facilitate the analysis of inhomogeneous systems, such as the one used in the present study, rather than homogeneous systems such as pure water. Because the number of trajectories of interest obtained by a single MD run is very limited, longer production runs are required for higher-accuracy in the case of simulations for inhomogeneous systems. Hence, the computational cost becomes a critical factor. Third, this full adaptive QM/MM method also provides access to large systems that cannot be treated with a full QM method. In the present study, we considered one hydronium ion and 2047 water molecules in a cubic box with side lengths of 39.5 Å, which is beyond the range of applications of ordinary AIMD simulations. It is known that the diffusion coefficient is subject to a limited size effect under periodic boundary condi-25 tions, resulting in a drastic underestimation in previous studies using AIMD simulations.</p><p>On the other hand, the SCMP method can suppress the artifacts, as we have previously reported, 41 although it is not sufficient. Therefore, in the present study, a major part of the deviation of the obtained diffusion coefficient from the experimental value can be attributed to either/both the shortcoming of the employed QM model (DFTB3/3OB) or/and the missing nuclear quantum effect. Owing to these advantages, the present full adaptive QM/MM method makes the hydronium ion simulation accessible with a plausible cost and moderate computation time, which has been a long-standing challenge in molecular simulations. It is clear that this will become an effective tool in advancing the theoretical analysis of hydronium ions in the subsequent stage.</p>