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 beneﬁt of accuracy found in the QM method and that of cost eﬃciency found in the MM method. However, it is diﬃcult to directly apply the QM/MM method to the dynamics of solution systems, particu-larly 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-ﬂy revisions on molecular deﬁnitions, is required for both the solute and solvent. Although several solvent-adaptive methods have been proposed, a full adaptive framework, which 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 conﬁrm that this method can stably and accurately simulate proton transfer dynamics in bulk water.


Introduction
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.
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-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. 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 2 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.
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).
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 QM/MM border, leading to the collapse of the MD simulation. Therefore, both definitions 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.
Potential energy and forces are independently evaluated for respective QM/MM partitions.
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: 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 where F (n) i 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.
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 where r i is the coordinate of the ith oxygen atom. 44 The weight function for the ith oxygen atom can be written as where r ij is the distance between the ith oxygen and the jth hydrogen atoms, and r o and d 5 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) 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 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 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.
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.

Theory and Method
Excess proton indicator In the present study, we propose a new EPI ξ, which was used as the QM center in the SCMP simulation, 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 8 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.
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.
Then, its Gaussian-type score function φ (r) is defined as where α and r 0 are the parameters (see the Discussion for details). Next, we assumed a normalized score φ iα as follows: Then, the virtual site η i for the ith oxygen atom was defined using the hydrogen coordinates r λ as 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: Here, S 1 satisfies the following boundary conditions: where R small and R large are the bonding parameters that satisfy R small < R large . Here, S 1 (r iα ) 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: Based on S 2 , the weight W i for the ith oxygen atom is defined as 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).
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.
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.
Here, the EPI ξ was adopted as the atom in the QM center. O where λ ( j (n) is the progress function of the respective jth QM or MM solvent, which continuously ranges between zero and one. I 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 where a, b ∈ {x, y, z}. If ξ({r}) is a function of the other variables in the Hamiltonian, then which does not satisfy condition 19. However, when ξ is an independent variable, Eq. 20 will be 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.
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 aṡ The MD algorithm for the ith particle is as follows: 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.
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.
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 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: 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: where δ 13 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.

Computational details
The SCMP method was implemented in a local version of the GROMACS 5.0.7 package. [51][52][53] In

Results
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,  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.
Indeed, MD simulations with larger values of c tend to become stuck earlier.
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  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.
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 18 that distance range.
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 δ.
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  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.
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 21 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.
However, it is difficult to precisely predict that it can indirectly have competing effects on water diffusion with respect to surrounding water molecules. 69

Discussions
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. Figure 8 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%.
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.
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.
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 23 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

Conclusions
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 func-