Search

We consider a range of model potentials with metastable states undergoing molecular dynamics coupled to a thermal bath in the high friction regime, and consider how the optimal reaction coordinate depends on the diffusion anisotropy. For this we use our recently proposed method 'Spectral gap optimization of order parameters (SGOOP)' (Tiwary and Berne, Proc. Natl. Acad. Sci. 113 2839 2016). We show how available information about dynamical observables in addition to static information can be incorporated into SGOOP, which can then be used to accurately determine the 'best' reaction coordinate for arbitrary anisotropies. We compare our results with transmission coefficient calculations and published benchmarks where applicable or available respectively.

We use a recently proposed method called Spectral Gap Optimization of Order Parameters (SGOOP) (Tiwary and Berne, Proc. Natl. Acad. Sci 2016, 113, 2839 (2016)), to determine an optimal 1-dimensional reaction coordinate (RC) for the unbinding of a bucky-ball from a pocket in explicit water. This RC is estimated as a linear combination of the multiple available order parameters that collectively can be used to distinguish the various stable states relevant for unbinding. We pay special attention to determining and quantifying the degree to which water molecules should be included in the RC. Using SGOOP with under-sampled biased simulations, we predict that water plays a distinct role in the reaction coordinate for unbinding in the case when the ligand is sterically constrained to move along an axis of symmetry. This prediction is validated through extensive calculations of the unbinding times through metadynamics, and by comparison through detailed balance with unbiased molecular dynamics estimate of the binding time. However when the steric constraint is removed, we find that the role of water in the reaction coordinate diminishes. Here instead SGOOP identifies a good one-dimensional RC involving various motional degrees of freedom., Comment: 7 pages, 5 figures

We consider the rate of transition for a particle between two metastable states coupled to a thermal environment for various magnitudes of the coupling strength, using the recently proposed infrequent metadynamics approach (Tiwary and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)). We are interested in understanding how this approach for obtaining rate constants performs as the dynamics regime changes from energy diffusion to spatial diffusion. Reassuringly, we find that the approach works remarkably well for various coupling strengths in the strong coupling regime, and to some extent even in the weak coupling regime., Comment: 3 pages, 1 figure, submitted to J. Chem. Phys

In modern day simulations of many-body systems much of the computational complexity is shifted to the identification of slowly changing molecular order parameters called collective variables (CV) or reaction coordinates. A vast array of enhanced sampling methods are based on the identification and biasing of these low-dimensional order parameters, whose fluctuations are important in driving rare events of interest. Here describe a new algorithm for finding optimal low-dimensional collective variables for use in enhanced sampling biasing methods like umbrella sampling, metadynamics and related methods, when limited prior static and dynamic information is known about the system, and a much larger set of candidate CVs is specified. The algorithm involves estimating the best combination of these candidate CVs, as quantified by a maximum path entropy estimate of the spectral gap for dynamics viewed as a function of that CV. Through multiple practical examples, we show how this post-processing procedure can lead to optimization of CV and several orders of magnitude improvement in the convergence of the free energy calculated through metadynamics, essentially giving the ability to extract useful information even from unsuccessful metadynamics runs., Comment: 7 pages, 4 figures; corrected missing figure number and added a reference

A key factor influencing a drug's efficacy is its residence time in the binding pocket of the host protein. Using atomistic computer simulation to predict this residence time and the associated dissociation process is a desirable but extremely difficult task due to the long timescales involved. This gets further complicated by the presence of biophysical factors such as steric and solvation effects. In this work, we perform molecular dynamics (MD) simulations of the unbinding of a popular prototypical hydrophobic cavity-ligand system using a metadynamics based approach that allows direct assessment of kinetic pathways and parameters. When constrained to move in an axial manner, we find the unbinding time to be on the order of 4000 sec. In accordance with previous studies, we find that the ligand must pass through a region of sharp dewetting transition manifested by sudden and high fluctuations in solvent density in the cavity. When we remove the steric constraints on ligand, the unbinding happens predominantly by an alternate pathway, where the unbinding becomes 20 times faster, and the sharp dewetting transition instead becomes continuous. We validate the unbinding timescales from metadynamics through a Poisson analysis, and by comparison through detailed balance to binding timescale estimates from unbiased MD. This work demonstrates that enhanced sampling can be used to perform explicit solvent molecular dynamics studies at timescales previously unattainable, obtaining direct and reliable pictures of the underlying physio-chemical factors including free energies and rate constants., Comment: 7 pages, 4 figures, supplementary PDF file, submitted

Heterogeneous ice growth exhibits a maximum in freezing rate arising from the competition between kinetics and the thermodynamic driving force between the solid and liquid states. Here, we use molecular dynamics simulations to elucidate the atomistic details of this competition, focusing on water properties in the interfacial region along the secondary prismatic direction. The crystal growth velocity is maximized when the efficiency of converting interfacial water molecules to ice, collectively known as the attachment kinetics, is greatest. We find water molecules that contact the intermediate ice layer in concave regions along the atomistically roughened surface are more likely to freeze directly. The increased roughening of the solid surface at large undercoolings consequently plays an important limiting role on the rate of ice growth, as water molecules are unable to integrate into increasingly deeper surface pockets. These results provide insights into the molecular mechanisms for self-assembly of solid phases that are important in many biological and atmospheric processes., Comment: Revised version 8 pages, 5 figures

When two phases of water are at equilibrium, the ratio of hydrogen isotopes in each is slightly altered due to their different phase affinities. This isotopic fractionation process can be utilized to analyze water's movement in the world's climate. Here we show that equilibrium fractionation ratios, an entirely quantum mechanical property, also provide a sensitive probe to assess the magnitude of nuclear quantum fluctuations in water. By comparing the predictions of a series of water models, we show that those describing the OH chemical bond as rigid or harmonic greatly over-predict the magnitude of isotope fractionation. Models that account for anharmonicity in this coordinate are shown to provide much more accurate results due to their ability to give partial cancellation between inter and intra-molecular quantum effects. These results give evidence of the existence of competing quantum effects in water and allow us to identify how this cancellation varies across a wide range of temperatures. In addition, this work demonstrates that simulation can provide accurate predictions and insights into hydrogen fractionation., Comment: 6 pages, 2 figures, 2 tables

Longstanding mechanistic questions about the role of protecting osmolyte trimethylamine N- oxide (TMAO) which favors protein folding and the denaturing osmolyte urea are addressed by studying their effects on the folding of uncharged polymer chains. Using atomistic molecular dynamics simulations, we show that 1-M TMAO and 7-M urea solutions act dramatically differently on these model polymer chains. Their behaviors are sensitive to the strength of the attractive dispersion interactions of the chain with its environment: when these dispersion interactions are high enough, TMAO suppresses the formation of extended conformations of the hydrophobic polymer as compared to water, while urea promotes formation of extended conformations. Similar trends are observed experimentally on real protein systems. Quite surprisingly, we find that both protecting and denaturing osmolytes strongly interact with the polymer, seemingly in contrast with existing explanations of the osmolyte effect on proteins. We show that what really matters for a protective osmolyte is its effective depletion as the polymer conformation changes, which leads to a negative change in the preferential binding coefficient. For TMAO, there is a much more favorable free energy of insertion of a single osmolyte near collapsed conformations of the polymer than near extended conformations. By contrast, urea is preferentially stabilized next to the extended conformation and thus has a denaturing effect.

A model of protein-ligand binding kinetics in which slow solvent dynamics results from hydrophobic drying transitions is investigated. Molecular dynamics simulations show that solvent in the receptor pocket can fluctuate between wet and dry states with lifetimes in each state that are long enough for the extraction of a separable potential of mean force and wet-to-dry transitions. We introduce a Diffusive Surface Hopping Model that is represented by a two-dimensional Markovian master equation. One dimension is the standard reaction coordinate, the ligand-pocket separation, and the other is the solvent state in the region between ligand and binding pocket which specifies whether it is wet or dry. In our model, the ligand diffuses on a dynamic free energy surface which undergoes kinetic transitions between the wet and dry states. The model yields good agreement with results from explicit solvent molecular dynamics simulation and an improved description of the kinetics of hydrophobic assembly. Furthermore, it is consistent with a "non-Markovian Brownian theory" for the ligand-pocket separation coordinate alone., Comment: 17 pages, 10 Figures

We study the kinetics of assembly of two plates of varying hydrophobicity, including cases where drying occurs and water strongly solvates the plate surfaces. The potential of mean force and molecular-scale hydrodynamics are computed from molecular dynamics simulations in explicit solvent as a function of particle separation. In agreement with our recent work on nanospheres [J. Phys. Chem. B 116, 378 (2012)] regions of high friction are found to be engendered by large and slow solvent fluctuations. These slow fluctuations can be due to either drying or confinement. The mean first passage times for assembly are computed by means of molecular dynamics simulations in explicit solvent and by Brownian dynamics simulations along the reaction path. Brownian dynamics makes use of the potential of mean force and hydrodynamic profile that we determined. Surprisingly, we find reasonable agreement between full scale molecular dynamics and Brownian dynamics, despite the role of slow solvent relaxation in the assembly process. We found that molecular scale hydrodynamic interactions are essential in describing the kinetics of assembly., Comment: 6 figures, 13 pages

Solvent plays an important role in the relative motion of nanoscopic bodies, and the study of such phenomena can help elucidate the mechanism of hydrophobic assembly, as well as the influence of solvent-mediated effects on in vivo motion in crowded cellular environments. Here we study important aspects of this problem within the framework of Brownian dynamics. We compute the free energy surface that the Brownian particles experience and their hydrodynamic interactions from molecular dynamics simulations in explicit solvent. We find that molecular scale effects dominate at short distances, thus giving rise to deviations from the predictions of continuum hydrodynamic theory. Drying phenomena, solvent layering, and fluctuations engender distinct signatures of the molecular scale. The rate of assembly in the diffusion-controlled limit is found to decrease from molecular scale hydrodynamic interactions, in opposition to the free energy driving force for hydrophobic assembly, and act to reinforce the influence of the free energy surface on the association of more hydrophilic bodies., Comment: 20 pages, 10 figures

The very nature of glass is somewhat mysterious: while relaxation times in glasses are of sufficient magnitude that large-scale motion on the atomic level is essentially as slow as it is in the crystalline state, the structure of glass appears barely different than that of the liquid that produced it. Quantum mechanical systems ranging from electron liquids to superfluid helium appear to form glasses, but as yet no unifying framework exists connecting classical and quantum regimes of vitrification. Here we develop new insights from theory and simulation into the quantum glass transition that surprisingly reveal distinct regions where quantum fluctuations can either promote or inhibit glass formation., Comment: Accepted for publication in Nature Physics. 22 pages, 3 figures, 1 Table

Multiple time scale molecular dynamics enhances computational efficiency by updating slow motions less frequently than fast motions. However, in practice the largest outer time step possible is limited not by the physical forces but by resonances between the fast and slow modes. In this paper we show that this problem can be alleviated by using a simple colored noise thermostatting scheme which selectively targets the high frequency modes in the system. For two sample problems, flexible water and solvated alanine dipeptide, we demonstrate that this allows the use of large outer time steps while still obtaining accurate sampling and minimizing the perturbation of the dynamics. Furthermore, this approach is shown to be comparable to constraining fast motions, thus providing an alternative to molecular dynamics with constraints., Comment: accepted for publication by the Journal of Chemical Physics

We study by molecular dynamics simulations the wetting/dewetting transition and the dependence of the free energy on distance between plates that contain both hydrophobic and hydrophilic particles. We show that dewetting and strength of hydrophobic interaction is very sensitive to the distribution of hydrophobic and hydrophilic domains. In particular, we find that plates characterized by a large domain of hydrophobic sites induce a dewetting transition and an attractive solvent-induced interaction. On the other hand, a homogeneous distribution of the hydrophobic and hydrophilic particles on the plates prevents the dewetting transition and produces a repulsive solvent-induced interaction. We also present results for a kind of Janus interface in which one plate consists of hydrophobic particles and the other of hydrophilic particles showing that the inter-plate gap remains wet until steric constraints at small separations eject the water molecules. Our results indicate that the Cassie equation, for the contact angle of a heterogeneous plate, can not be used to predict the critical distance of dewetting. These results indicate that hydrophobic interactions between nanoscale surfaces with strong large length-scale hydrophobicity can be highly cooperative and thus they argue against additivity of the hydrophobic interactions between different surface domains in these cases. These findings are pertinent to certain protein-protein interactions where additivity is commonly assumed., Comment: 28 pages, 6 figures

We propose a general methodology for calculating the self-diffusion tensor from molecular dynamics for a liquid with a liquid-gas or liquid-solid interface. The standard method used in bulk fluids, based on computing the mean square displacement as a function of time and extracting the asymptotic linear time dependence from this, is not valid for systems with interfaces or for confined fluids. The method proposed here is based on imposing virtual boundary conditions on the molecular system and computing survival probabilities and specified time correlation functions in different layers of the fluid up to and including the interfacial layer. By running dual simulations, one based on MD and the other based on Langevin dynamics, using the same boundary conditions, one can fit the Langevin survival probability at long times to the MD computed survival probability, thereby determining the diffusion coefficient as a function of distance of the layers from the interface. We compute the elements of the diffusion tensor of water as a function of distance from the liquid vapor interface of water. Far from the interface the diffusion tensor is found to be isotropic, as expected, and the diffusion coefficient has the value $D\approx$ .22\AA$^2$/psec in agreement with what is found in the bulk liquid. In the interfacial region the diffusion tensor is axially anisotropic, with values of $D_{\parallel}\approx$. 8\AA$^2$/psec and $D_{\perp}\approx$. 5\AA$^2$/psec for the components parallel and normal the interface surface respectively. We also show that diffusion in confined geometries can be calculated by imposing appropriate boundary conditions on the molecular system and computing time correlation functions of the eigenfunctions of the diffusion operator corresponding to the same boundary conditions.

A new molecular dynamics model in which the point charges on atomic sites are allowed to fluctuate in response to the environment is developed and applied to water. The idea for treating charges as variables is based on the concept of electronegativity equalization according to which: (a) The electronegativity of an atomic site is dependent on the atom's type and charge and is perturbed by the electrostatic potential it experiences from its neighbors and (b) Charge is transferred between atomic sites in such a way that electronegativities are equalized. The charges are treated as dynamical variables using an extended Lagrangian method in which the charges are given a fictitious mass, velocities and kinetic energy and then propagated according to Newtonian mechanics along with the atomic degrees of freedom. Models for water with fluctuating charges are developed using the geometries of two common fixed-charge water potentials: the simple point charge (SPC) and the 4-point transferable intermolecular potential (TIP4P). Both fluctuating charge models give accurate predictions for gas-phase and liquid state properties, including radial distribution functions, the dielectric constant, and the diffusion constant. The method does not introduce any new intermolecular interactions beyond those already present in the fixed charge models and increases the computer time by only a factor of 1.1, making this method tractable for large systems., Comment: J. Chem. Phys., in press 29 pages + 12 figures

The calculation of the solvation properties of a single water molecule in liquid water is carried out in two ways. In the first, the water molecule is placed in a cavity and the solvent is treated as a dielectric continuum. This model is analyzed by numerically solving the Poisson equation using the DelPhi program. The resulting solvation properties depend sensitively on the shape and size of the cavity. In the second method, the solvent and solute molecules are treated explicitly in molecular dynamics simulations using Ewald boundary conditions. We find a 2 kcal/mole difference in solvation free energies predicted by these two methods when standard cavity radii are used. In addition, dielectric continuum theory assumes that the solvent reacts solely by realigning its electric moments linearly with the strength of the solute's electric field; the results of the molecular simulation show important non-linear effects. Non-linear solvent effects are generally of two types: dielectric saturation, due to solvent-solute hydrogen bonds, and electrostriction, a decrease in the solute cavity due to an increased electrostatic interaction. We find very good agreement between the two methods if the radii defining the solute cavity used in the continuum theory is decreased with the solute charges, Comment: J. of the American Chemical Society