A polarizable multistate empirical valence bond model for proton transport in aqueous solution

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.

An effective partial charge model for bulk and surface properties of cubic ZrO2, Y2O3 and yttrium-stabilised zirconia

In this work a newly parametrised Coulomb plus Buckingham potential formulation for cubic ZrO₂, Y₂O₃ and yttrium-stabilised zirconia (YSZ) is presented. The density and pair distributions obtained for neat ZrO₂ and Y₂O₃ under ambient conditions are in excellent agreement with experimental data, while the vibrational power spectra are highly similar compared to those obtained via ab initio molecular dynamics simulations at the PBEsol level. In addition, it is shown that the use of effective partial charges has several advantages compared to interaction potentials employing the oxidation states in the evaluation of the coulombic interactions: (i) the diffusion coefficient and the associated activation energy of oxygen ions evaluated for YSZn (n = 4 to 12) display the best agreement with experimental data; (ii) no unphysical reorganisation of the interface and the bulk are observed in simulations of the (110) and (111) surfaces of cubic ZrO₂ and Y₂O₃, while due to the strong coulombic contributions in the case of the tested full-charge models a pronounced restructuring of the interface and the bulk is observed in the ZrO₂ case, and (iii) the use of effective partial charges ensures compatibility with existing solvent models and force-fields for the treatment of molecular compounds.

Ion Pairing in HCl-Water Clusters: From Electronic Structure Investigations to Multiconfigurational Force-Field Development

In the bulk, condensed-phase HCl exists as a dissociated Cl^- ion and a proton that is delocalized over solvating water molecules. However, in the gas phase, HCl is covalent, and even on the introduction of hydrating water molecules, the HCl covalent state dominates small clusters and is relevant at larger clusters including 21 water molecules. Electronic structure calculations (at the MP2 level) and ab initio metadynamics simulations (at the DFT level) have been carried out on HCl-(H₂O)_n clusters with n = 2-22 to investigate distinct solvation environments in clusters from covalent HCl structure, to contact ion pairs and solvent-separated ion pairs. The data were further used to train and validate a multiconfigurational force-field for HCl-water clusters that incorporates covalent HCl states into the MS-EVB3.2 formalism. Additionally, the many-body interaction of the Cl^- ion with water and the excess proton was modeled by the introduction of two geometric three-body terms that incorporates the dominant many-body interaction in an efficient noniterative manner.

Maxwell-Boltzmann versus non-ergodic events in the velocity distribution of water molecules evaporated from protonated water nanodroplets

Measurement of velocity distributions of evaporated water monomers from small mass- and energy-selected protonated water clusters allows probing the extent of thermalization after excitation of these ultimately small nanodroplets. Electronic excitation of a molecule in the cluster is here induced by a single collision with an argon atom in the keV energy range. The measured velocity distributions of the departing neutral molecules exhibit bimodal shapes with a lower-velocity part consistent with a complete redistribution of the deposited energy in the entire cluster and a higher-velocity contribution corresponding to evaporation before complete energy redistribution. Statistical molecular dynamics calculations reproduce the bimodal shape of the velocity distributions by assuming an initial spreading of the excitation energy among all modes, thereby reproducing the lower velocity contribution of the distribution. By contrast, assuming the deposited energy to be initially localized among the modes of a single molecule leads to calculated distributions with two components whose shape is in accordance with the experimental results. The characteristics and the relative abundance of these two contributions in the velocity distributions obtained are presented and discussed as a function of the number of molecules (n = 2-10) in the ionized nanodroplet H⁺(H₂O) _n .

Towards a dissociative SPC-like water model II. The impact of Lennard-Jones and Buckingham non-coulombic forces

The dissociative water potential by Garofalini and coworkers has been re-formulated in the framework of the widely employed Lennard-Jones and Buckingham potentials, enhancing the transferability of the model to third party simulation programs.

A comparison of sodium and hydrogen halides at the air-water interface

New molecular models, parameterized to ab initio calculations, were developed to describe HBr and HI at the air-water interface. These were used to compare how the air-water interface influenced dissociation of NaX and HX, with X being Cl, Br, or I, and also their propensity for the interface. The polarizable multistate empirical valence bond method, which explicitly describes proton sharing, was used to model HX. Results showed that the air-water interface suppressed HX dissociation from a contact ion pair to a solvent separated to a greater degree than NaX dissociation. Furthermore, HX had a greater propensity for the interface than NaX, which was a consequence of the hydronium ion having a greatest interfacial activity of all species studied. As a consequence of this, the average configuration of dissociated HX, while in both contact ion and solvent separated ion pairs near the air-water interface, is with the dissociated hydrogen oriented more towards the air than the X atom.

Restrained Proton Indicator in Combined Quantum-Mechanics/Molecular-Mechanics Dynamics Simulations of Proton Transfer through a Carbon Nanotube

Recently, a collective variable "proton indicator" was purposed for tracking an excess proton solvated in bulk water in molecular dynamics simulations. In this work, we demonstrate the feasibility of utilizing the position of this proton indicator as a reaction coordinate to model an excess proton migrating through a hydrophobic carbon nanotube in combined quantum-mechanics/molecular-mechanics simulations. Our results indicate that applying a harmonic restraint to the proton indicator in the bulk solvent near the nanotube pore entrance leads to the recruitment of water molecules into the pore. This is consistent with an earlier study that employed a multistate empirical valence bond potential and a different representation (center of excess charge) of the proton. We attribute this water recruitment to the delocalized nature of the solvated proton, which prefers to be in high-dielectric bulk solvent. While water recruitment into the pore is considered an artifact in the present simulations (because of the artificially imposed restraint on the proton), if the proton were naturally restrained, it could assist in building water wires prior to proton transfer through the pore. The potential of mean force for a proton translocation through the water-filled pore was computed by umbrella sampling, where the bias potentials were applied to the proton indicator. The free energy curve and barrier heights agree reasonably with those in the literature. The results suggest that the proton indicator can be used as a reaction coordinate in simulations of proton transport in confined environments.

Potential energy construction in the diabatic picture for quantum mechanical rate constants of intermolecular proton transfer

We propose a simple method for potential construction in the diabatic picture and the estimation of thermal rate constants for intermolecular proton transfer reactions using quantum dynamics simulations carried out on the constructed potentials.

Towards a dissociative SPC-like water model – probing the impact of intramolecular Coulombic contributions

A modification of the dissociative Garofalini water model towards an SPC-like Coulombic formulation proved to enhance accuracy and transferability of this successful force field approach.

Velocity of a Molecule Evaporated from a Water Nanodroplet: Maxwell-Boltzmann Statistics versus Non-Ergodic Events

The velocity of a molecule evaporated from a mass-selected protonated water nanodroplet is measured by velocity map imaging in combination with a recently developed mass spectrometry technique. The measured velocity distributions allow probing statistical energy redistribution in ultimately small water nanodroplets after ultrafast electronic excitation. As the droplet size increases, the velocity distribution rapidly approaches the behavior expected for macroscopic droplets. However, a distinct high-velocity contribution provides evidence of molecular evaporation before complete energy redistribution, corresponding to non-ergodic events.

Role of Presolvation and Anharmonicity in Aqueous Phase Hydrated Proton Solvation and Transport

Results from condensed phase ab initio molecular dynamics (AIMD) simulations suggest a proton transfer reaction is facilitated by "presolvation" in which the hydronium is transiently solvated by four water molecules, similar to the typical solvation structure of water, by accepting a weak hydrogen bond from the fourth water molecule. A new version 3.2 multistate empirical valence bond (MS-EVB 3.2) model for the hydrated excess proton incorporating this presolvation behavior is therefore developed. The classical MS-EVB simulations show similar structural properties as those of the previous model but with significantly improved diffusive behavior. The inclusion of nuclear quantum effects in the MS-EVB also provides an even better description of the proton diffusion rate. To quantify the influence of anharmonicity, a second model (aMS-EVB 3.2) is developed using the anharmonic aSPC/Fw water model, which provides similar structural properties but improved spectroscopic responses at high frequencies.

Tuning the reactivity of a dissociative force field: proton transfer properties of aqueous H3O+ and their dependence on the three-body interaction

Selective adjustment of the three-body interaction of a dissociative water potential results in a significant improvement in the description of proton transport properties.

Double layer effects in a model of proton discharge on charged electrodes

We report first results on double layer effects on proton discharge reactions from aqueous solutions to charged platinum electrodes. We have extended a recently developed combined proton transfer/proton discharge model on the basis of empirical valence bond theory to include specifically adsorbed sodium cations and chloride anions. For each of four studied systems 800-1000 trajectories of a discharging proton were integrated by molecular dynamics simulations until discharge occurred. The results show significant influences of ion presence on the average behavior of protons prior to the discharge event. Rationalization of the observed behavior cannot be based solely on the electrochemical potential (or surface charge) but needs to resort to the molecular details of the double layer structure.

Perspectives for hybrid ab initio/molecular mechanical simulations of solutions: from complex chemistry to proton-transfer reactions and interfaces

As a consequence of the ongoing development of enhanced computational resources, theoretical chemistry has become an increasingly valuable field for the investigation of a variety of chemical systems. Simulations employing a hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) technique have been shown to be a particularly promising approach, whenever ultrafast (i.e., picosecond) dynamical properties are to be studied, which are in many cases difficult to access via experimental techniques. Details of the quantum mechanical charge field (QMCF) ansatz, an advanced QM/MM protocol, are discussed and simulation results for various systems ranging from simple ionic hydrates to solvated organic molecules and coordination complexes in solution are presented. A particularly challenging application is the description of proton-transfer reactions in chemical simulations, which is a prerequisite to study acidified and basic systems. The methodical requirements for a combination of the QMCF methodology with a dissociative potential model for the description of the solvent are discussed. Furthermore, the possible extension of QM/MM approaches to solid/liquid interfaces is outlined.

Role of the Membrane Dipole Potential for Proton Transport in Gramicidin A Embedded in a DMPC Bilayer

The membrane potential at the water/phospholipid interfaces may play a key role for proton conduction of gramicidin A (gA). Here we address this issue by Density Functional Theory-based molecular dynamics and metadynamics simulations. The calculations, performed on gA embedded in a solvated 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) model membrane environment (about 2,000 atoms), indicate that (i) the membrane dipole potential rises at the channel mouth by ∼0.4 V. A similar value has been measured for gA embedded in a DMPC monolayer; (ii) the calculated free energy barrier is located at the channel entrance, consistent with experiments comparing gA proton conduction in different bilayers. The electronic structures of the proton ligands (water molecules and peptide units) are similar to those in the bulk solvent. Based on these results, we suggest an important role of the membrane dipole potential for the free energy barrier of proton permeation of gA. This may provide a rationale for the large increase in the rate of proton conduction under application of a transmembrane voltage, as observed experimentally. Our calculations might suggest also a role for proton desolvation for the permeation process. This role has already emerged from EVB calculations on gA embedded in a model membrane.

Lifetimes of excess protons in water using a dissociative water potential

Molecular dynamics simulations using a dissociative water potential were applied to study transport of excess protons in water and determine the applicability of this potential to describe such behavior. While originally developed for gas-phase molecules and bulk liquid water, the potential is transferrable to nanoconfinement and interface scenarios. Applied here, it shows proton behavior consistent with ab initio calculations and empirical models specifically designed to describe proton transport. Both Eigen and Zundel complexes are observed in the simulations showing the Eigen-Zundel-Eigen-type mechanism. In addition to reproducing the short-time rattling of the excess proton between the two oxygens of Zundel complexes, a picosecond-scale lifetime was also found. These longer-lived H3O(+) ions are caused by the rapid conversion of the local solvation structure around the transferring proton from a Zundel-like form to an Eigen-like form following the transfer, effectively severing the path along which the proton can rattle. The migration of H(+) over long times (>100 ps) deviates from the conventional short-time multiexponentially decaying lifetime autocorrelation model and follows the t(-3/2) power-law behavior. The potential function employed here matches many of the features of proton transport observed in ab initio molecular dynamics simulations as well as the highly developed empirical valence bond models, yet is computationally very efficient, enabling longer time and larger systems to be studied.

Multiscale reactive molecular dynamics

Many processes important to chemistry, materials science, and biology cannot be described without considering electronic and nuclear-level dynamics and their coupling to slower, cooperative motions of the system. These inherently multiscale problems require computationally efficient and accurate methods to converge statistical properties. In this paper, a method is presented that uses data directly from condensed phase ab initio simulations to develop reactive molecular dynamics models that do not require predefined empirical functions. Instead, the interactions used in the reactive model are expressed as linear combinations of interpolating functions that are optimized by using a linear least-squares algorithm. One notable benefit of the procedure outlined here is the capability to minimize the number of parameters requiring nonlinear optimization. The method presented can be generally applied to multiscale problems and is demonstrated by generating reactive models for the hydrated excess proton and hydroxide ion based directly on condensed phase ab initio molecular dynamics simulations. The resulting models faithfully reproduce the water-ion structural properties and diffusion constants from the ab initio simulations. Additionally, the free energy profiles for proton transfer, which is sensitive to the structural diffusion of both ions in water, are reproduced. The high fidelity of these models to ab initio simulations will permit accurate modeling of general chemical reactions in condensed phase systems with computational efficiency orders of magnitudes greater than currently possible with ab initio simulation methods, thus facilitating a proper statistical sampling of the coupling to slow, large-scale motions of the system.

Natural polarizability and flexibility via explicit valency: the case of water

As the dominant physiological solvent, water drives the folding of biological macromolecules, influences conformational changes, determines the ionization states of surface groups, actively participates in catalytic events, and provides "wires" for long-range proton transfer. Elucidation of all these roles calls for atomistic simulations. However, currently available methods do not lend themselves to efficient simulation of proton transfer events, or even polarizability and flexibility. Here, we report that an explicit account of valency can provide a unified description for the polarizability, flexibility, and dissociability of water in one intuitive and efficient setting. We call this approach LEWIS, after the chemical theory that inspires the use of valence electron pairs. In this paper, we provide details of the method, the choice of the training set, and predictions for the neat ambient liquid, with emphasis on structure, dynamics, and polarization. LEWIS water provides a good description of bulk properties, and dipolar and quadrupolar responses.

Correlation effects on the structure and dynamics of the H3O+ hydrate: B3LYP/MM and MP2/MM MD simulations

Two combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations, namely B3LYP/MM and MP2/MM, have been performed to investigate the possible influence of electron correlation on the structure and dynamics of the H(3)O(+) hydrate. In comparison to the previously published HF/MM results, both B3LYP/MM and MP2/MM simulations clearly reveal stronger H(3)O(+)-water hydrogen bond interactions, which are reflected in a slightly greater compactness of the H(3)O(+) hydrate. However, the B3LYP/MM simulation, although providing structural details very close to the MP2/MM data, shows an artificially slow dynamic nature of some first shell water molecules as a consequence of the formation of a long-lived H(3)O(+)···H(2)O hydrogen bonding structure.

Constant pH Molecular Dynamics in Explicit Solvent with λ-Dynamics

pH is an important parameter in condensed-phase systems, because it determines the protonation state of titratable groups and thus influences the structure, dynamics, and function of molecules in solution. In most force field simulation protocols, however, the protonation state of a system (rather than its pH) is kept fixed and cannot adapt to changes of the local environment. Here, we present a method, implemented within the MD package GROMACS, for constant pH molecular dynamics simulations in explicit solvent that is based on the λ-dynamics approach. In the latter, the dynamics of the titration coordinate λ, which interpolates between the protonated and deprotonated states, is driven by generalized forces between the protonated and deprotonated states. The hydration free energy, as a function of pH, is included to facilitate constant pH simulations. The protonation states of titratable groups are allowed to change dynamically during a simulation, thus reproducing average protonation probabilities at a certain pH. The accuracy of the method is tested against titration curves of single amino acids and a dipeptide in explicit solvent.

Hydration structure in dilute hydrofluoric acid

We have performed the multistate empirical valence bond (MS-EVB) molecular dynamics simulations of a dilute hydrofluoric acid solution at ambient temperature to study the hydration structure associated with its weak acidity. The developed MS-EVB model showed reasonable agreement with experimental and previous ab initio molecular dynamics and reference interaction site model self-consistent field simulations for the free energy and structural properties. The local tetrahedral and translational order parameters around the fluorine atom significantly increase in the transition and product states of the HF dissociation reaction. This indicates that the angular and translational rearrangements of the hydrogen-bond topology are necessary especially around the fluorine atom. At the transition state of the proton transfer, the tetrahedral order parameters are very large, whereas the translational order parameters are not. This suggests that for the proton transfer to occur the large angular rearrangements of the hydrogen-bond topology are more necessary than the translational ones.

Proton transfer to charged platinum electrodes. A molecular dynamics trajectory study

A recently developed empirical valence bond (EVB) model for proton transfer on Pt(111) electrodes (Wilhelm et al 2008 J. Phys. Chem. C 112 10814) has been applied in molecular dynamics (MD) simulations of a water film in contact with a charged Pt surface. A total of seven negative surface charge densities σ between -7.5 and -18.9 µC cm(-2) were investigated. For each value of σ, between 30 and 84 initial conditions of a solvated proton within a water slab were sampled, and the trajectories were integrated until discharge of a proton occurred on the charged surfaces. We have calculated the mean rates for discharge and for adsorption of solvated protons within the adsorbed water layer in contact with the metal electrode as a function of surface charge density. For the less negative values of σ we observe a Tafel-like exponential increase of discharge rate with decreasing σ. At the more negative values this exponential increase levels off and the discharge process is apparently transport limited. Mechanistically, the Tafel regime corresponds to a stepwise proton transfer: first, a proton is transferred from the bulk into the contact water layer, which is followed by transfer of a proton to the charged surface and concomitant discharge. At the more negative surface charge densities the proton transfer into the contact water layer and the transfer of another proton to the surface and its discharge occur almost simultaneously.

A multistate empirical valence bond model for solvation and transport simulations of OH- in aqueous solutions

We describe a new multistate empirical valence bond (MS-EVB) model of OH(-) in aqueous solutions. This model is based on the recently proposed "charged ring" parameterization for the intermolecular interaction of hydroxyl ion with water [Ufimtsev, et al., Chem. Phys. Lett., 2007, 442, 128] and is suitable for classical molecular simulations of OH(-) solvation and transport. The model reproduces the hydration structure of OH(-)(aq) in good agreement with experimental data and the results of ab initio molecular dynamics simulations. It also accurately captures the major structural, energetic, and dynamic aspects of the proton transfer processes involving OH(-) (aq). The model predicts an approximately two-fold increase of the OH(-) mobility due to proton exchange reactions.

Infrared spectra of a model phenol-amine proton transfer complex in nanoconfined CH3Cl

The vibrational spectra of a model phenol-amine proton transfer complex dissolved in CH3Cl solvent confined in a 12 A radius spherical hydrophobic cavity were calculated using mixed quantum-classical molecular dynamics simulations. The reaction free energy of the proton transfer complex was varied in order to explore the contributions to the vibrational absorption band from product and reactant species. The vibrational spectra of the model proton transfer complex resulted in motionally narrowed spectral linewidths with two distinct peaks for products and reactants in cases where the system undergoes chemical exchange. It was found that the n=1 and n=2 vibrational excited states combine to form diabatic states such that the spectra have contributions from both n=0 --> n=1 and n=0 --> n=2 transitions. A strong relationship between the instantaneous vibrational frequency and a collective solvent coordinate was found that assists in understanding the origin of the spectral features.

Acid-base analysis: a critique of the Stewart and bicarbonate-centered approaches

When approaching the analysis of disorders of acid-base balance, physical chemists, physiologists, and clinicians, tend to focus on different aspects of the relevant phenomenology. The physical chemist focuses on a quantitative understanding of proton hydration and aqueous proton transfer reactions that alter the acidity of a given solution. The physiologist focuses on molecular, cellular, and whole organ transport processes that modulate the acidity of a given body fluid compartment. The clinician emphasizes the diagnosis, clinical causes, and most appropriate treatment of acid-base disturbances. Historically, two different conceptual frameworks have evolved among clinicians and physiologists for interpreting acid-base phenomena. The traditional or bicarbonate-centered framework relies quantitatively on the Henderson-Hasselbalch equation, whereas the Stewart or strong ion approach utilizes either the original Stewart equation or its simplified version derived by Constable. In this review, the concepts underlying the bicarbonate-centered and Stewart formulations are analyzed in detail, emphasizing the differences in how each approach characterizes acid-base phenomenology at the molecular level, tissue level, and in the clinical realm. A quantitative comparison of the equations that are currently used in the literature to calculate H+concentration ([H+]) is included to clear up some of the misconceptions that currently exist in this area. Our analysis demonstrates that while the principle of electroneutrality plays a central role in the strong ion formulation, electroneutrality mechanistically does not dictate a specific [H+], and the strong ion and bicarbonate-centered approaches are quantitatively identical even in the presence of nonbicarbonate buffers. Finally, our analysis indicates that the bicarbonate-centered approach utilizing the Henderson-Hasselbalch equation is a mechanistic formulation that reflects the underlying acid-base phenomenology.

Molecular Dynamics Analysis of Interfacial Structures and Sum Frequency Generation Spectra of Aqueous Hydrogen Halide Solutions

Molecular dynamics simulation for gas/liquid interfaces of aqueous hydrochloric (HCl) and hydroiodic (HI) acid solutions is performed to calculate and analyze their sum frequency generation (SFG) spectra. The present MD simulation supports the strong preference of hydronium ions at the topmost surface layer and a consequent formation of ionic double layers by the hydronium and halide ions near the interface. Accordingly, the orientational order of surface water in the double layers is reversed in the acid solutions from that in the salt (NaCl or NaI). The calculated SFG spectra of the O-H stretching region reproduce the experimental spectra of ssp and sps polarizations well. In the ssp spectra, the strong enhancement in the hydrogen-bonding region for the acid solutions is elucidated by two mechanisms, ordered orientation of water in the double layer and symmetric OH stretching of the surface hydronium ions. In the sps spectra, reversed orientation of surface water is evidenced in the spectral line shapes, which are quite different from those of the salt solutions.

Ab initio QM/MM dynamics of H3O+ in water

A molecular dynamics (MD) simulation based on a combined ab initio quantum mechanics/molecular mechanics (QM/MM) method has been performed to investigate the solvation structure and dynamics of H3O+ in water. The QM region is a sphere around the central H3O+ ion, and contains about 6-8 water molecules. It is treated at the Hartree-Fock (HF) level, while the rest of the system is described by means of classical pair potentials. The Eigen complex (H9O4+) is found to be the most prevalent species in the aqueous solution, partly due to the selection scheme of the center of the QM region. The QM/MM results show that the Eigen complex frequently converts back and forth into the Zundel (H5O2+) structure. Besides the three nearest-neighbor water molecules directly hydrogen-bonded to H3O+, other neighbor waters, such as a fourth water molecule which interacts preferentially with the oxygen atom of the hydronium ion, are found occasionally near the ion. Analyses of the water exchange processes and the mean residence times of water molecules in the ion's hydration shell indicate that such next-nearest neighbor water molecules participate in the rearrangement of the hydrogen bond network during fluctuative formation of the Zundel ion and, thus, contribute to the Grotthuss transport of the proton.