How Short Is Too Short for the Interactions of a Water Potential? Exploring the Parameter Space of a Coarse-Grained Water Model Using Uncertainty Quantification

Accelerated Bayesian Inference for Molecular Simulations using Local Gaussian Process Surrogate Models

While Bayesian inference is the gold standard for uncertainty quantification and propagation, its use within physical chemistry encounters formidable computational barriers. These bottlenecks are magnified for modeling data with many independent variables, such as X-ray/neutron scattering patterns and electromagnetic spectra. To address this challenge, we employ local Gaussian process (LGP) surrogate models to accelerate Bayesian optimization over these complex thermophysical properties. The time-complexity of the LGPs scales linearly in the number of independent variables, in stark contrast to the computationally expensive cubic scaling of conventional Gaussian processes. To illustrate the method, we trained a LGP surrogate model on the radial distribution function of liquid neon and observed a 1,760,000-fold speed-up compared to molecular dynamics simulation, beating a conventional GP by three orders-of-magnitude. We conclude that LGPs are robust and efficient surrogate models poised to expand the application of Bayesian inference in molecular simulations to a broad spectrum of experimental data.

Reverse Mapping of Coarse Grained Polyglutamine Conformations from PRIME20 Sampling

An inverse coarse-graining protocol is presented for generating and validating atomistic structures of large (bio-) molecules from conformations obtained via a coarse-grained sampling method. Specifically, the protocol is implemented and tested based on the (coarse-grained) PRIME20 protein model (P20/SAMC), and the resulting all-atom conformations are simulated using conventional biomolecular force fields. The phase space sampling at the coarse-grained level is performed with a stochastical approximation Monte Carlo approach. The method is applied to a series of polypeptides, specifically dimers of polyglutamine with varying chain length in aqueous solution. The majority (>70 %) of the conformations obtained from the coarse-grained peptide model can successfully be mapped back to atomistic structures that remain conformationally stable during 10 ns of molecular dynamics simulations. This work can be seen as the first step towards the overarching goal of improving our understanding of protein aggregation phenomena through simulation methods.

Kinetics and Mechanisms of Pressure-Induced Ice Amorphization and Polyamorphic Transitions in a Machine-Learned Coarse-Grained Water Model

Water glasses have attracted considerable attention due to their potential connection to a liquid-liquid transition in supercooled water. Here we use molecular simulations to investigate the formation and phase behavior of water glasses using the machine-learned bond-order parameter (ML-BOP) water model. We produce glasses through hyperquenching of water, pressure-induced amorphization (PIA) of ice, and pressure-induced polyamorphic transformations. We find that PIA of polycrystalline ice occurs at a lower pressure than that of monocrystalline ice and through a different mechanism. The temperature dependence of the amorphization pressure of polycrystalline ice for ML-BOP agrees with that in experiments. We also find that ML-BOP accurately reproduces the density, coordination number, and structural features of low-density (LDA), high-density (HDA), and very high-density (VHDA) amorphous water glasses. ML-BOP accurately reproduces the experimental radial distribution function of LDA but overpredicts the minimum between the first two shells in high-density glasses. We examine the kinetics and mechanism of the transformation between low-density and high-density glasses and find that the sharp nature of these transitions in ML-BOP is similar to that in experiments and all-atom water models with a liquid-liquid transition. Transitions between ML-BOP glasses occur through a spinodal-like mechanism, similar to ice crystallization from LDA. Both glass-to-glass and glass-to-ice transformations have Avrami-Kolmogorov kinetics with exponent n = 1.5 ± 0.2 in experiments and simulations. Importantly, ML-BOP reproduces the competition between crystallization and HDA→LDA transition above the glass transition temperature T_g, and separation of their time scales below T_g, observed also in experiments. These findings demonstrate the ability of ML-BOP to accurately reproduce water properties across various regimes, making it a promising model for addressing the competition between polyamorphic transitions and crystallization in water and solutions.

Stability and Metastability of Liquid Water in a Machine-Learned Coarse-Grained Model with Short-Range Interactions

Coarse-grained water models are ∼100 times more efficient than all-atom models, enabling simulations of supercooled water and crystallization. The machine-learned monatomic model ML-BOP reproduces the experimental equation of state (EOS) and ice-liquid thermodynamics at 0.1 MPa on par with the all-atom TIP4P/2005 and TIP4P/Ice models. These all-atom models were parametrized using high-pressure experimental data and are either accurate for water's EOS (TIP4P/2005) or ice-liquid equilibrium (TIP4P/Ice). ML-BOP was parametrized from temperature-dependent ice and liquid experimental densities and melting data at 0.1 MPa; its only pressure training is from compression of TIP4P/2005 ice at 0 K. Here we investigate whether ML-BOP replicates the experimental EOS and ice-water thermodynamics along all pressures of ice I. We find that ML-BOP reproduces the temperature, enthalpy, entropy, and volume of melting of hexagonal ice up to 400 MPa and the EOS of water along the melting line with an accuracy that rivals that of both TIP4P/2005 and TIP4P/Ice. We interpret that the accuracy of ML-BOP originates from its ability to capture the shift between compact and open local structures to changes in pressure and temperature. ML-BOP reproduces the sharpening of the tetrahedral peak of the pair distribution function of water upon supercooling, and its pressure dependence. We characterize the region of metastability of liquid ML-BOP with respect to crystallization and cavitation. The accessibility of ice crystallization to simulations of ML-BOP, together with its accurate representation of the thermodynamics of water, makes it promising for investigating the interplay between anomalies, glass transition, and crystallization under conditions challenging to access through experiments.

Oblique Impacts of Nanodroplets upon Surfaces

In this work, oblique impacts of nanodroplets impacting surfaces in a wide range of impact angles (α) are investigated in detail via molecular dynamics simulations. Five outcomes are observed, including deposition, prompt splashing, break-up, separation, and shattering. With increasing impact angle, the outcomes of prompt splashing, break-up, separation, and shattering are enlarged but the one of deposition is compressed. By drawing a We_n ∼ α phase diagram, the outcome regimes and corresponding boundaries of them can be successfully identified, and the boundary between the deposition and other outcome regimes is theoretically modeled and shows good agreement with the phase diagram, where We_n is the normal impact Weber number. For further understanding of the oblique impacts, the maximum spreading factor, as the feature parameter of spreading, is investigated. Asymmetry spreading behaviors are observed, noting that β_max,∥ is always larger than β_max,⊥. β_max,⊥ is tested that it only depends on We_n with wide impact angles and could be predicted by the scaling law of β_max,⊥ = 0.7We_n^1/4. However, β_max,∥ depends on not only We_n but also impact angles. A modified model is proposed for predicting β_max,∥ as 0.7We_n^1/4 + 0.001(We_n tan² α)^3/2, which shows good agreement with data on surfaces with θ from 73 to 105° in wide We_n and α ranges.

Is Ice Nucleation by Organic Crystals Nonclassical? An Assessment of the Monolayer Hypothesis of Ice Nucleation

Potent ice nucleating organic crystals display an increase in nucleation efficiency with pressure and memory effect after pressurization that set them apart from inorganic nucleants. These characteristics were proposed to arise from an ordered water monolayer at the organic-water interface. It was interpreted that ordering of the monolayer is the limiting step for ice nucleation on organic crystals, rendering their mechanism of nucleation nonclassical. Despite the importance of organics in atmospheric ice nucleation, that explanation has never been investigated. Here we elucidate the structure of interfacial water and its role in ice nucleation at ambient pressure on phloroglucinol dihydrate, the paradigmatic example of outstanding ice nucleating organic crystal, using molecular simulations. The simulations confirm the existence of an interfacial monolayer that orders on cooling and becomes fully ordered upon ice formation. The monolayer does not resemble any ice face but seamlessly connects the distinct hydrogen-bonding orders of ice and the organic surface. Although large ordered patches develop in the monolayer before ice nucleates, we find that the critical step is the formation of the ice crystallite, indicating that the mechanism is classical. We predict that the fully ordered, crystalline monolayer nucleates ice above -2 °C and could be responsible for the exceptional ice nucleation by the organic crystal at high pressures. The lifetime of the fully ordered monolayer around 0 °C, however, is too short to account for the memory effect reported in the experiments. The latter could arise from an increase in the melting temperature of ice confined by strongly ice-binding surfaces.

Uncertainty Quantification and Sensitivity Analysis of Partial Charges on Macroscopic Solvent Properties in Molecular Dynamics Simulations with a Machine Learning Model

The molecular dynamics (MD) simulation technique is among the most broadly used computational methods to investigate atomistic phenomena in a variety of chemical and biological systems. One of the most common (and most uncertain) parametrization steps in MD simulations of soft materials is the assignment of partial charges to atoms. Here, we apply uncertainty quantification and sensitivity analysis calculations to assess the uncertainty associated with partial charge assignment in the context of MD simulations of an organic solvent. Our results indicate that the effect of partial charge variance on bulk properties, such as solubility parameters, diffusivity, dipole moment, and density, measured from MD simulations is significant; however, measured properties are observed to be less sensitive to partial charges of less accessible (or buried) atoms. Diffusivity, for example, exhibits a global sensitivity of up to 22 × 10^-5 cm²/s per electron charge on some acetonitrile atoms. We then demonstrate that machine learning techniques, such as Gaussian process regression (GPR), can be effective and rapid tools for uncertainty quantification of MD simulations. We show that the formulation and application of an efficient GPR surrogate model for the prediction of responses effectively reduces the computational time of additional sample points from hours to milliseconds. This study provides a much-needed context for the effect that partial charge uncertainty has on MD-derived material properties to illustrate the benefit of considering partial charges as distributions rather than point-values. To aid in this treatment, this work then demonstrates methods for rapid characterization of resulting sensitivity in MD simulations.

Model reduction of rigid-body molecular dynamics via generalized multipole potentials

Motivated by the challenges of uncertainty quantification for coarse-grained (CG) molecular dynamics, we investigate the role of perturbation theory in model reduction of classical systems. In particular, we consider the task of coarse-graining rigid bodies in the context of generalized multipole potentials that have controllable levels of accuracy relative to their atomistic counterparts. We show how the multipole framework yields a hierarchy of models that systematically connects a CG "point molecule" approximation to the exact dynamics. We use these results to understand when and how the CG models fail to describe atomistic dynamics at the trajectory level and develop asymptotic error estimates for approximate molecular potential energies. Implications for other model-reduction strategies are also discussed. Key findings of this work are that (i) omitting rotational energy introduces significant error when coarse-graining and (ii) attention to symmetry can improve accuracy of "point-molecule" approximations. Analytical derivations and numerical results support these conclusions. Relevance to nonrigid bodies is also discussed.

How Do Surfactants Control the Agglomeration of Clathrate Hydrates?

Clathrate hydrates can spontaneously form under typical conditions found in oil and gas pipelines. The agglomeration of clathrates into large solid masses plugs the pipelines, posing adverse safety, economic, and environmental threats. Surfactants are customarily used to prevent the aggregation of clathrate particles and their coalescence with water droplets. It is generally assumed that a large contact angle between the surfactant-covered clathrate and water is a key predictor of the antiagglomerant performance of the surfactant. Here we use molecular dynamic simulations to investigate the structure and dynamics of surfactant films at the clathrate-oil interface, and their impact on the contact angle and coalescence between water droplets and hydrate particles. In agreement with the experiments, the simulations predict that surfactant-covered clathrate-oil interfaces are oil wet but super-hydrophobic to water. Although the water contact angle determines the driving force for coalescence, we find that a large contact angle is not sufficient to predict good antiagglomerant performance of a surfactant. We conclude that the length of the surfactant molecules, the density of the interfacial film, and the strength of binding of its molecules to the clathrate surface are the main factors in preventing the coalescence and agglomeration of clathrate particles with water droplets in oil. Our analysis provides a molecular foundation to guide the molecular design of effective clathrate antiagglomerants.

Machine learning coarse grained models for water

An accurate and computationally efficient molecular level description of mesoscopic behavior of ice-water systems remains a major challenge. Here, we introduce a set of machine-learned coarse-grained (CG) models (ML-BOP, ML-BOP_dih, and ML-mW) that accurately describe the structure and thermodynamic anomalies of both water and ice at mesoscopic scales, all at two orders of magnitude cheaper computational cost than existing atomistic models. In a significant departure from conventional force-field fitting, we use a multilevel evolutionary strategy that trains CG models against not just energetics from first-principles and experiments but also temperature-dependent properties inferred from on-the-fly molecular dynamics (~ 10's of milliseconds of overall trajectories). Our ML BOP models predict both the correct experimental melting point of ice and the temperature of maximum density of liquid water that remained elusive to-date. Our ML workflow navigates efficiently through the high-dimensional parameter space to even improve upon existing high-quality CG models (e.g. mW model).

Bayesian selection for coarse-grained models of liquid water

The necessity for accurate and computationally efficient representations of water in atomistic simulations that can span biologically relevant timescales has born the necessity of coarse-grained (CG) modeling. Despite numerous advances, CG water models rely mostly on a-priori specified assumptions. How these assumptions affect the model accuracy, efficiency, and in particular transferability, has not been systematically investigated. Here we propose a data driven comparison and selection for CG water models through a Hierarchical Bayesian framework. We examine CG water models that differ in their level of coarse-graining, structure, and number of interaction sites. We find that the importance of electrostatic interactions for the physical system under consideration is a dominant criterion for the model selection. Multi-site models are favored, unless the effects of water in electrostatic screening are not relevant, in which case the single site model is preferred due to its computational savings. The charge distribution is found to play an important role in the multi-site model’s accuracy while the flexibility of the bonds/angles may only slightly improve the models. Furthermore, we find significant variations in the computational cost of these models. We present a data informed rationale for the selection of CG water models and provide guidance for future water model designs.

What Controls the Limit of Supercooling and Superheating of Pinned Ice Surfaces?

Cold-adapted organisms produce antifreeze proteins and glycoproteins to control the growth, melting and recrystallization of ice. It has been proposed that these molecules pin the crystal surface, creating a curvature that arrests the growth and melting of the crystal. Here we use thermodynamic modeling and molecular simulations to demonstrate that the curvature of the superheated or supercooled surface depends on the temperature and distances between ice-binding molecules, but not the details of their interactions with ice. We perform simulations of ice pinned with the antifreeze protein TmAFP, polyvinyl alcohol with different degrees of polymerization, and model ice-binding molecules to determine the thermal hystereses on melting and freezing, i.e. the maximum curvature that can be attained before, respectively, ice melts or grows irreversibly over the ice-binding molecules. We find that the thermal hysteresis is controlled by the bulkiness of the ice-binding molecules and their footprint at the ice surface. We elucidate the origin of the asymmetry between freezing and melting hysteresis found in experiments and propose guidelines to design synthetic antifreeze molecules with potent thermal hysteresis activity.

Ice-Nucleating and Antifreeze Proteins Recognize Ice through a Diversity of Anchored Clathrate and Ice-like Motifs

Cold-adapted organisms produce antifreeze and ice-nucleating proteins to prevent and promote ice formation. The crystal structure of hyperactive bacterial antifreeze protein (AFP) MpAFP suggests that this protein binds ice through an anchored clathrate motif. It is not known whether other hyperactive AFPs and ice-nucleating proteins (INPs) use the same motif to recognize or nucleate ice. Here we use molecular simulations to elucidate the ice-binding motifs of hyperactive insect AFPs and a model INP of Pseudomonas syringae. We find that insect AFPs recognize ice through anchored clathrate motifs distinct from that of MpAFP. By performing simulations of ice nucleation by PsINP, we identify two distinct ice-binding sites on opposite sides of the β-helix. The ice-nucleating sequences identified in the simulations agree with those previously proposed for the closely related INP of Pseudomonas borealis based on the structure of the protein. The simulations indicate that these sites have comparable ice-nucleating efficiency, but distinct binding motifs, controlled by the amino acid sequence: one is an anchored clathrate and the other ice-like. We conclude that anchored clathrate and ice-like motifs can be equally effective for binding proteins to ice and promoting ice nucleation.

Is Water at the Graphite Interface Vapor-like or Ice-like?

Graphitic surfaces are the main component of soot, a major constituent of atmospheric aerosols. Experiments indicate that soots of different origins display a wide range of abilities to heterogeneously nucleate ice. The ability of pure graphite to nucleate ice in experiments, however, seems to be almost negligible. Nevertheless, molecular simulations with the monatomic water model mW with water-carbon interactions parameterized to reproduce the experimental contact angle of water on graphite predict that pure graphite nucleates ice. According to classical nucleation theory, the ability of a surface to nucleate ice is controlled by the binding free energy between ice immersed in liquid water and the surface. To establish whether the discrepancy in freezing efficiencies of graphite in mW simulations and experiments arises from the coarse resolution of the model or can be fixed by reparameterization, it is important to elucidate the contributions of the water-graphite, water-ice, and ice-water interfaces to the free energy, enthalpy, and entropy of binding for both water and the model. Here we use thermodynamic analysis and free energy calculations to determine these interfacial properties. We demonstrate that liquid water at the graphite interface is not ice-like or vapor-like: it has similar free energy, entropy, and enthalpy as water in the bulk. The thermodynamics of the water-graphite interface is well reproduced by the mW model. We find that the entropy of binding between graphite and ice is positive and dominated, in both experiments and simulations, by the favorable entropy of reducing the ice-water interface. Our analysis indicates that the discrepancy in freezing efficiencies of graphite in experiments and the simulations with mW arises from the inability of the model to simultaneously reproduce the contact angle of liquid water on graphite and the free energy of the ice-graphite interface. This transferability issue is intrinsic to the resolution of the model, and arises from its lack of rotational degrees of freedom.

Bayesian calibration of coarse-grained forces: Efficiently addressing transferability

Generating and calibrating forces that are transferable across a range of state-points remains a challenging task in coarse-grained (CG) molecular dynamics. In this work, we present a coarse-graining workflow, inspired by ideas from uncertainty quantification and numerical analysis, to address this problem. The key idea behind our approach is to introduce a Bayesian correction algorithm that uses functional derivatives of CG simulations to rapidly and inexpensively recalibrate initial estimates f0 of forces anchored by standard methods such as force-matching. Taking density-temperature relationships as a running example, we demonstrate that this algorithm, in concert with various interpolation schemes, can be used to efficiently compute physically reasonable force curves on a fine grid of state-points. Importantly, we show that our workflow is robust to several choices available to the modeler, including the interpolation schemes and tools used to construct f0. In a related vein, we also demonstrate that our approach can speed up coarse-graining by reducing the number of atomistic simulations needed as inputs to standard methods for generating CG forces.

Parameterization of a coarse-grained model with short-ranged interactions for modeling fuel cell membranes with controlled water uptake

The coarse-grained model FF_pvap reproduces the experimental activity coefficient of water in tetramethylammonium chloride solutions over a wide range of concentrations, with a hundred-fold gain in computing efficiency with respect to atomistic models.

Flexible Boundaries for Multiresolution Solvation: An Algorithm for Spatial Multiscaling in Molecular Dynamics Simulations

An algorithm is proposed for performing molecular dynamics (MD) simulations of a biomolecular solute represented at atomistic resolution surrounded by a surface layer of atomistic fine-grained (FG) solvent molecules within a bulk represented by coarse-grained (CG) solvent beads. The method, called flexible boundaries for multiresolution solvation (FBMS), is based on: (i) a three-region layering of the solvent around the solute, involving an FG layer surrounded by a mixed FG-CG buffer layer, itself surrounded by a bulk CG region; (ii) a definition of the layer boundary that relies on an effective distance to the solute surface and is thus adapted to the shape of the solute as well as adjusts to its conformational changes. The effective surface distance is defined by inverse-nth power averaging over the distances to all non-hydrogen solute atoms, and the layering is enforced by means of half-harmonic distance restraints, attractive for the FG molecules and repulsive for the CG beads. A restraint-free region at intermediate distances enables the formation of the buffer layer, where the FG and CG solvents can mix freely. The algorithm is tested and validated using the GROMOS force field and the associated FG (SPC) and CG (polarizable CGW) water models. The test systems include pure-water systems, where one FG molecule plays the role of a solute, and a deca-alanine peptide with two widely different solute shapes considered, α-helical and fully extended. In particular, as the peptide unfolds, the number of FG molecules required to fill its close-range solvation layer increases, with the additional molecules being provided by the buffer layer. Further validation involves simulations of four proteins in multiresolution FG/CG mixtures. The resulting structural, energetic, and solvation properties are found to be similar to those observed in corresponding pure FG simulations.

Triplet correlation functions in liquid water

Triplet correlations have been shown to play a crucial role in the transformation of simple liquids to anomalous tetrahedral fluids [M. Singh, D. Dhabal, A. H. Nguyen, V. Molinero, and C. Chakravarty, Phys. Rev. Lett. 112, 147801 (2014)]. Here we examine triplet correlation functions for water, arguably the most important tetrahedral liquid, under ambient conditions, using configurational ensembles derived from molecular dynamics (MD) simulations and reverse Monte Carlo (RMC) datasets fitted to experimental scattering data. Four different RMC data sets with widely varying hydrogen-bond topologies fitted to neutron and x-ray scattering data are considered [K. T. Wikfeldt, M. Leetmaa, M. P. Ljungberg, A. Nilsson, and L. G. M. Pettersson, J. Phys. Chem. B 113, 6246 (2009)]. Molecular dynamics simulations are performed for two rigid-body effective pair potentials (SPC/E and TIP4P/2005) and the monatomic water (mW) model. Triplet correlation functions are compared with other structural measures for tetrahedrality, such as the O-O-O angular distribution function and the local tetrahedral order distributions. In contrast to the pair correlation functions, which are identical for all the RMC ensembles, the O-O-O triplet correlation function can discriminate between ensembles with different degrees of tetrahedral network formation with the maximally symmetric, tetrahedral SYM dataset displaying distinct signatures of tetrahedrality similar to those obtained from atomistic simulations of the SPC/E model. Triplet correlations from the RMC datasets conform closely to the Kirkwood superposition approximation, while those from MD simulations show deviations within the first two neighbour shells. The possibilities for experimental estimation of triplet correlations of water and other tetrahedral liquids are discussed.

Coarse-Graining of TIP4P/2005, TIP4P-Ew, SPC/E, and TIP3P to Monatomic Anisotropic Water Models Using Relative Entropy Minimization

Coarse-grained models are becoming a competitive alternative for modeling processes that occur over time and length scales beyond the reach of fully atomistic molecular simulations. Ideally, coarse-grained models should not only achieve high computational efficiency but also provide accurate predictions and fundamental insight into the role of molecular interactions, the characteristic behavior, and properties of the system they model. In this work we derive a series of monatomic coarse-grained water models mX(REM) from the most popular atomistic water models X = TIP3P, SPC/E, TIP4P-Ew, and TIP4P/2005, using the relative entropy minimization (REM) method. Each coarse-grained water molecule is represented by a single particle that interacts through short-ranged anisotropic interactions that encourage the formation of "hydrogen-bonded" structures. We systematically investigate the features of the coarse-grained models in reproducing over 20 structural, dynamic, and thermodynamic properties of the reference atomistic water models-including the existence and locus of the characteristic density anomaly. The mX(REM) coarse-grained models reproduce quite faithfully the radial and angular distribution function of water, produce a temperature of maximum density (TMD), and stabilize the ice I crystal. Moreover, the ratio between the TMD and the melting temperature of the crystal in the mX(REM) models and liquid-ice equilibrium properties show reasonable agreement with the results of the corresponding atomistic models. The mX(REM) models, however, severely underestimate the cohesive energy of the condensed water phases. We investigate which specific limitations of the coarse-grained models arise from the REM methodology, from the monatomic nature of the models, and from the Stillinger-Weber interaction potential form. Our analysis indicates that a small compromise in the accuracy of structural properties can result in a significant increase of the overall accuracy and representability of the coarse-grained water models. We evaluate the accuracy of the atomistic and the monatomic anisotropic coarse-grained water models, including the mW water model, in reproducing experimental water properties. We find that mW and mTIP4P/2005(REM) score closer to experiment than widely used atomistic water models. We conclude that monatomic models of water with short-range, anisotropic "hydrogen-bonding" three-body interactions can be competitive in accuracy with fully atomistic models for the study of a wide range of properties and phenomena at less than 1/100th of the computational cost.