Biasing the Center of Charge in Molecular Dynamics Simulations with Empirical Valence Bond Models:  Free Energetics of an Excess Proton in a Water Droplet

Abstract

Multistate empirical valence bond (EVB) models provide an accurate description of the energetics of proton transfer and solvation in complex molecular systems and can be efficiently used in molecular dynamics computer simulations. Within such models, the location of the moving protonic charge can be specified by the so-called center of charge, defined as a weighted average over the diabatic states of the EVB model. In this paper, we use first-order perturbation theory to calculate the molecular forces that arise if a bias potential is applied to the center of charge. Such bias potentials are often necessary when molecular dynamics simulations are used to determine free energies related to proton transfer and not all relevant proton positions are sampled with sufficient frequency during the available computing time. The force expressions we derive are easy to evaluate and do not create any significant computational cost compared with unbiased EVB simulations. As an illustration of the method, we study proton transfer in a small liquid water droplet consisting of 128 water molecules plus an excess proton. Contrary to predictions of continuum electrostatics, but in agreement with previous computer simulations of similar systems, we observe that the excess proton is predominantly located at the surface of the droplet. Using the formalism developed in this paper, we calculate the reversible work required to carry the protonic charge from the droplet surface to its core, finding a value of roughly 4 k(B)T.

Water at hydrophobic interfaces delays proton surface-to-bulk transfer and provides a pathway for lateral proton diffusion

“…gests that a free energy barrier ΔG H þ is intrinsic to any interface separating hydrophobic and hydrophilic phases. This supports the opinions from previous simulations (15)(16)(17)(18)(19)(20)(21). Proton retention close to the surface decreases interfacial pH at the first two water layers by about 2-3 units below bulk pH, as indicated by our extensive 75 ps-long ab initio metadynamics-based free energy simulations, including 1,707 atoms.…”

“…5) because, beyond that, G bind ¼ 0 (Fig. 7A); (iii) since our calculated value of the free energy barrier ΔG H þ ¼ 6 AE 2 RT is not very different from those calculated for water/ vapor (17,18), water/carbon nanotube (19), and water∕CCl 4 (20) interfaces, some of these findings may be general to dielectric mismatched interfaces.…”

“…G Born-image . This contribution arises from moving a proton from the hydrophilic environment of bulk water (ε approximately 80) to the more hydrophobic environment adjacent to the interface (ε approximately [10][11][12][13][14][15][16][17][18][19][20]. We use continuum electrostatics for the calculation of this contribution (SI Text).…”

“…For example, for the water/vapor interface, ab initio and MS-EVB simulations revealed the acidification of the top surface water layer as compared to bulk water (15)(16)(17)(18). Free energy minima about 3 Å wide from the interface were also found in MS-EVB simulations at water/carbon nanotubes (19) and water∕ CCl 4 interfaces (20).…”

“…According to simulations carried out with nonbiological models (15)(16)(17)(18)(19)(20)(21), well-defined binding sites may be not required for proton attraction to interfaces between high and low dielectric regions. For example, for the water/vapor interface, ab initio and MS-EVB simulations revealed the acidification of the top surface water layer as compared to bulk water (15)(16)(17)(18).…”