Improved Indicator Algorithms for Tracking a Hydrated Proton as A Local Structural Defect in Grotthuss Diffusion in Aqueous Solutions

Multi-proton dynamics near membrane-water interface

Protons are crucial for biological energy transduction between membrane proteins. While experiments suggest rapid proton motion over large distances at the membrane-water interface, computational studies employing a single excess proton found the proton immobilized near the lipid headgroup. To address this discrepancy, we conduct DFTB3 simulations by incrementally adding protons up to three. We show that a single proton moves rapidly toward the nearest headgroup, where it is either repelled by a choline group or binds covalently to phosphatic oxygen. With multiple protons, while some are trapped by the lipid headgroups, the remaining proton diffuses laterally faster than in bulk water. Driven by excess energy, this proton initially jumps to the center of the water slab before relaxing into the third- and second-hydration shells. Lateral diffusion rates increase as the proton stabilizes in the second hydration shell. These results provide insights into proton dynamics near membranes and explain experimental observations.

Reshaping the QM Region On-the-Fly: Adaptive-Shape QM/MM Dynamic Simulations of a Hydrated Proton in Bulk Water

Adaptive quantum mechanics/molecular mechanics (QM/MM) reclassifies on-the-fly a molecule or molecular fragment as QM or MM during dynamics simulations without abrupt changes in the energy or forces. Notably, the permuted adaptive-partitioning (PAP) algorithms have been applied to simulate a hydrated proton, with a mobile QM zone anchored at a pseudoatom called a proton indicator. The position of the proton indicator approximates the location of the delocalized excess proton, yielding a smooth trajectory of the proton diffusing via the Grotthuss mechanism in aqueous solutions. The mobile QM zone, which has been taken to be a sphere with a preset radius, follows the proton wherever it goes. Although the simulations are successful, the use of a spherical QM zone has one disadvantage: A large preset radius must be utilized to minimize the chance of missing water molecules that are important to proton translocation. A large radius leads to a large QM zone, which is computationally expensive. In this work, we report a new way to set up the QM zone, where one includes only the water molecules important to proton transfer. The importance of a given water molecule is quantified by its "weight" that depends on its relation to the reaction path of proton transfer. The weight varies smoothly, ensuring that a water molecule gradually appears in or disappears from the QM zone without abrupt changes, as required by the PAP method. Consequently, the shape of the QM zone evolves on-the-fly, keeping the QM zone as small as possible and as large as necessary. Test simulations demonstrate that the new algorithm significantly improves the computation efficiency while maintaining the proper descriptions of proton transfer in bulk water.

Adaptive-Partitioning Multilayer Dynamics Simulations: 2. Implementations of the Permuted and Interpolated Adaptive-Partitioning Gradients

Recently, an adaptive-partitioning multilayer Q1/Q2/MM method was proposed, where Q1 and Q2 denote, respectively, two distinct quantum-mechanical levels of theory and MM, the molecular-mechanical force fields. Such a multilayer model resembles the ONIOM (our own N-layered integrated molecular orbital and molecular mechanics) model by Morokuma and co-workers, but it is distinguished by on-the-fly reclassifying atoms to be Q1, Q2, or MM in dynamics simulations. To smoothly blend the levels of descriptions of the atoms, buffer zones are introduced between adjacent layers, and the energy is smoothly interpolated. In particular, the Q1/Q2 interaction energy was expressed in two different formalisms: permuted and interpolated adaptive-partitioning (PAP and IAP), respectively. While the PAP energy is based on a weighted many-body expansion, the IAP energy is derived via alchemical quantum calculations with interpolated Fock and overlap matrices. In this article, we examine in-depth the irregularities in the IAP energy near the boundary between the buffer and Q2 zones, which were found prominent in some calculations. These irregularities are due to basis-set linear dependencies, which can be effectively suppressed using a cutoff for the weighted atomic orbital coefficients. Furthermore, we derived and implemented the gradients for both PAP and IAP. Test calculations on a series of water cluster models show perfectly smooth gradients in PAP, while a minor discontinuity occurs in IAP gradients at the buffer/Q2 boundary. The energy and gradient discontinuities in IAP become smaller when moving the buffer/Q2 boundary further away from the Q1 center and when increasing the size of the basis sets used. Overall, those discontinuities are controllable, and possible ways to further diminish them are discussed.

Tracking the Delocalized Proton in Concerted Proton Transfer in Bulk Water

A solvated proton in water is often characterized as a charge or structural defect, and it is important to track its evolution on-the-fly in certain dynamics simulations. Previously, we introduced the proton indicator, a pseudo-atom, whose position approximates the location of the excess proton modeled as a structural defect. The proton indicator generally yields a smooth trajectory of a hydrated proton diffusing in aqueous solutions, including in the events of stepwise proton transfer via the Grotthuss mechanism; however, the proton indicator did not perform well in the events of concerted proton transfer, for which it occasionally yielded large position displacements between two successive time steps. To overcome this hurdle, we develop a new algorithm of a proton indicator with greatly enhanced performance for concerted proton transfer in bulk water. A protocol is proposed to exhaustively explore the hydrogen-bonding network of the water wires over which the excess proton is delocalized and to properly account for the contributions of the water molecules in this network as the geometry evolves. The new proton indicator (called Indicator 2.0) is assessed in dynamics simulations of an excess proton in bulk water and in specially constructed model systems of more complex architectures. The results demonstrate that the new indicator yields a smooth trajectory in both stepwise and concerted proton transfers.

Grotthuss Molecular Dynamics Simulations for Modeling Proton Hopping in Electrosprayed Water Droplets

Excess protons in water exhibit unique transport properties because they can rapidly hop along H-bonded water wires. Considerable progress has been made in unraveling this Grotthuss diffusion mechanism using quantum mechanical-based computational techniques. Unfortunately, high computational cost tends to restrict those techniques to small systems and short times. Molecular dynamics (MD) simulations can be applied to much larger systems and longer time windows. However, standard MD methods do not permit the dissociation/formation of covalent bonds, such that Grotthuss diffusion cannot be captured. Here, we bridge this gap by combining atomistic MD simulations (using Gromacs and TIP4P/2005 water) with proton hopping. Excess protons are modeled as hydronium ions that undergo H₃O⁺ + H₂O → H₂O + H₃O⁺ transitions. In accordance with ab initio MD data, these Grotthuss hopping events are executed in "bursts" with quasi-instantaneous hopping across one or more waters. The bursts are separated by regular MD periods during which H₃O⁺ ions undergo Brownian diffusion. The resulting proton diffusion coefficient agrees with the literature value. We apply this Grotthuss MD technique to highly charged water droplets that are in a size regime encountered during electrospray ionization (5 nm radius, ∼17,000 H₂O). The droplets undergo rapid solvent evaporation and occasional H₃O⁺ ejection, keeping them at ca. 81% of the Rayleigh limit. The simulated behavior is consistent with phase Doppler anemometry data. The Grotthuss MD technique developed here should be useful for modeling the behavior of various proton-containing systems that are too large for high-level computational approaches. In particular, we envision future applications related to electrospray processes, where earlier simulations used metal cations while in reality excess protons dominate.