The control of assembly and crystallization of molecules is becoming increasingly important in chemistry, engineering, and materials sciences. Crystallization is also central to understand natural processes that include the formation of atmospheric ice and biomineralization. Organic surfaces, biomolecules, and even liquid/vapor interfaces can promote the nucleation of crystals. These soft surfaces present significant structural fluctuations, which have been shown to strongly impact the rate of crystallization. This raises the question of whether degrees of freedom of soft surfaces play a role in the reaction coordinate for crystal nucleation. Here we use molecular simulations to investigate the mechanism of ice nucleation promoted by an alcohol monolayer. Our analysis indicates that while the flexibility of the surface strongly depresses its ice nucleation ability, it does not play a role in the coordinate that controls the transformation from liquid to ice. We find that the variable that drives the transformation is the size of the crystalline cluster, the same as that for the homogeneous crystallization. We argue that this is a general result that arises from the separation of time scales between surface fluctuations and the crossing of the transition state barrier for crystallization.

Heterogeneous nucleation of ice induced by organic materials is of fundamental importance for climate, biology, and industry. Among organic ice-nucleating surfaces, monolayers of long chain alcohols are particularly effective, while monolayers of fatty acids are significantly less so. As these monolayers expose to water hydroxyl groups with an order that resembles the one in the basal plane of ice, it was proposed that lattice matching between ice and the surface controls their ice-nucleating efficiency. Organic monolayers are soft materials and display significant fluctuations. It has been conjectured that these fluctuations assist in the nucleation of ice. Here we use molecular dynamic simulations and laboratory experiments to investigate the relationship between the structure and fluctuations of hydroxylated organic surfaces and the temperature at which they nucleate ice. We find that these surfaces order interfacial water to form domains with ice-like order that are the birthplace of ice. Both mismatch and fluctuations decrease the size of the preordered domains and monotonously decrease the ice freezing temperature. The simulations indicate that fluctuations depress the freezing efficiency of monolayers of alcohols or acids to half the value predicted from lattice mismatch alone. The model captures the experimental trend in freezing efficiencies as a function of chain length and predicts that alcohols have higher freezing efficiency than acids of the same chain length. These trends are mostly controlled by the modulation of the structural mismatch to ice. We use classical nucleation theory to show that the freezing efficiencies of the monolayers are directly related to their free energy of binding to ice. This study provides a general framework to relate the equilibrium thermodynamics of ice binding to a surface and the nonequilibrium ice freezing temperature and suggests that these could be predicted from the structure of interfacial water.

The ice–air interface is an important locus of environmental chemical reactions. The structure and dynamics of the ice surface impact the uptake of trace gases and kinetics of reactions in the atmosphere and snowpack. At tropospheric temperatures, the ice surface is partially premelted. Experiments indicate that ions increase the liquidity of the ice surface but hydrophilic organics do not. However, it is not yet known the extent of the perturbation solutes induce at the ice surface and what is the role of the disordered liquid-like layer in modulating the interaction between solutes and their mobility and aggregation at the ice surface. Here we use large-scale molecular simulations to investigate the effect of ions and glyoxal, one of the most abundant oxygenated volatile organic compounds in the atmosphere, on the structure, dynamics, and solvation properties of the ice surface. We find that the premelted surface of ice has unique solvation properties, different from those of liquid water. The increase in surface liquidity resulting from the hydration of ions leads to a water-mediated attraction of ions at the ice surface. Glyoxal molecules, on the other hand, perturb only slightly the surface of ice and do not experience water-driven attraction. They nonetheless accumulate as dry agglomerates at the ice surface, driven by direct interactions between the organic molecules. The enhanced attraction and clustering of ions and organics at the ice surface may play a significant role in modulating the mechanism and rate of heterogeneous chemical reactions occurring at the surface of atmospheric ice particles.

The design of polymer electrolyte membranes with controlled water uptake is of high importance for high-performance fuel cells, because the water content of the membranes modulates their conductivity, chemical stability and mechanical strength. The water activity aw controls the equilibrium water uptake of a system. Predicting aw of materials is currently a daunting challenge for molecular simulations, because calculations of water activity require grand canonical simulations that are extremely expensive even with classical non-polarizable force fields. Moreover, force fields do not generally reproduce aw of solutions. Here, we first present a general strategy to parameterize force fields that reproduce the experimental aw of solutions, and then implement that strategy to re-parameterize the interactions in FFcomp, a coarse-grained model for hydrated polyphenylene oxide/trimethylamine chloride (PPO/TMACl) membranes in which the TMA cation is attached to the PPO backbone and the Cl anion is in the mobile water nanophase. Coarse-grained models based on short-ranged potentials successfully model fuel cell membranes and other concentrated aqueous electrolyte solutions because electrostatic interactions are highly screened in these systems. The new force field, FFpvap, differs from the original FFcomp only in the parameters of the ion–ion interactions, yet it reproduces aw in TMACl solutions with accuracy within 0.5 and 3% of the experimental value in all the concentration range relevant to the operation of fuel cell membranes. We find that the heat needed to vaporize water in solutions with as little as five water molecules per ion pair is essentially the same as in pure water, despite the strong water–ion interactions and their impact on the water activity. We review the literature to demonstrate that this is independent of the model and a general feature of water solutions. FFpvap reproduces the radial distribution functions and captures well the relative diffusivities of water and ions in the ionic solution as predicted by the reference atomistic GAFF-TIP4P/2005 model, while providing a hundred-fold gain in computing efficiency with respect to the atomistic model. With the backbone fragments inherited from FFcomp, the new FFpvap force field can be used to model hydrated polymer electrolyte membranes and advance the design of fuel cell membranes with controlled water uptake and conductivity.

Molecular simulations provide a versatile tool to study the structure, anion conductivity, and stability of anion-exchange membrane (AEM) materials and can provide a fundamental understanding of the relation between structure and property of membranes that is key for their use in fuel cells and other applications. The quest for large spatial and temporal scales required to model the multiscale structure and transport processes in the polymer electrolyte membranes, however, cannot be met with fully atomistic models, and the available coarse-grained (CG) models suffer from several challenges associated with their low-resolution. Here, we develop a high-resolution CG force field for hydrated polyphenylene oxide/trimethylamine chloride (PPO/TMACl) membranes compatible with the mW water model using a hierarchical parametrization approach based on Uncertainty Quantification and reference atomistic simulations modeled with the Generalized Amber Force Field (GAFF) and TIP4P/2005 water. The parametrization weighs multiple properties, including coordination numbers, radial distribution functions (RDFs), self-diffusion coefficients of water and ions, relative vapor pressure of water in the solution, hydration enthalpy of the tetramethylammonium chloride (TMACl) salt, and cohesive energy of its aqueous solutions. We analyze the interdependence between properties and address how to compromise between the accuracies of the properties to achieve an overall best representability. Our optimized CG model FFcomp quantitatively reproduces the diffusivities and RDFs of the reference atomistic model and qualitatively reproduces the experimental relative vapor pressure of water in solutions of tetramethylammonium chloride. These properties are of utmost relevance for the design and operation of fuel cell membranes. To our knowledge, this is the first CG model that includes explicitly each water and ion and accounts for hydrophobic, ionic, and intramolecular interactions explicitly parametrized to reproduce multiple properties of interest for hydrated polyelectrolyte membranes. The CG model of hydrated PPO/TMACl water is about 100 times faster than the reference atomistic GAFF-TIP4P/2005 model. The strategy implemented here can be used in the parametrization of CG models for other substances, such as biomolecular systems and membranes for desalination, water purification, and redox flow batteries. We anticipate that the large spatial and temporal simulations made possible by the CG model will advance the quest for anion-exchange membranes with improved transport and mechanical properties.

Ice crystals in the atmosphere nucleate from supercooled liquid water and grow by vapor uptake. The structure of the ice polymorph grown has strong impact on the morphology and light scattering of the ice crystals, modulates the amount of water vapor in ice clouds, and can impact the molecular uptake and reactivity of atmospheric aerosols. Experiments and molecular simulations indicate that ice nucleated and grown from deeply supercooled liquid water is metastable stacking disordered ice. The ice polymorph grown from vapor has not yet been determined. Here we use large-scale molecular simulations to determine the structure of ice that grows as a result of uptake of water vapor in the temperature range relevant to cirrus and mixed-phase clouds, elucidate the molecular mechanism of the formation of ice at the vapor interface, and compute the free energy difference between cubic and hexagonal ice interfaces with vapor. We find that vapor deposition results in growth of stacking disordered ice only under conditions of extreme supersaturation, for which a nonequilibrium liquid layer completely wets the surface of ice. Such extreme conditions have been used to produce stacking disordered frost ice in experiments and may be plausible in the summer polar mesosphere. Growth of ice from vapor at moderate supersaturations in the temperature range relevant to cirrus and mixed-phase clouds, from 200 to 260 K, produces exclusively the stable hexagonal ice polymorph. Cubic ice is disfavored with respect to hexagonal ice not only by a small penalty in the bulk free energy (3.6 ± 1.5 J mol–1 at 260 K) but also by a large free energy penalty at the ice–vapor interface (89.7 ± 12.8 J mol–1 at 260 K). The latter originates in higher vibrational entropy of the hexagonal-terminated ice–vapor interface. We predict that the free energy penalty against the cubic ice interface should decrease strongly with temperature, resulting in some degree of stacking disorder in ice grown from vapor in the tropical tropopause layer, and in polar stratospheric and noctilucent clouds. Our findings support and explain the evolution of the morphology of ice crystals from hexagonal to trigonal symmetry with decreasing temperature, as reported by experiments and in situ measurements in clouds. We conclude that selective growth of the elusive cubic ice polymorph by manipulation of the interfacial properties can likely be achieved at the ice–liquid interface but not at the ice–vapor interface.

The vapor pressure of water is a key property in a large class of applications from the design of membranes for fuel cells and separations to the prediction of the mixing state of atmospheric aerosols. Molecular simulations have been used to compute vapor pressures, and a few studies on liquid mixtures and solutions have been reported on the basis of the Gibbs Ensemble Monte Carlo method in combination with atomistic force fields. These simulations are costly, making them impractical for the prediction of the vapor pressure of complex materials. The goal of the present work is twofold: (1) to demonstrate the use of the grand canonical screening approach (Factorovich, M. H. J. Chem. Phys. 2014, 140, 064111) to compute the vapor pressure of solutions and to extend the methodology for the treatment of systems without a liquid–vapor interface and (2) to investigate the ability of computationally efficient high-resolution coarse-grained models based on the mW monatomic water potential and ions described exclusively with short-range interactions to reproduce the relative vapor pressure of aqueous solutions. We find that coarse-grained models of LiCl and NaCl solutions faithfully reproduce the experimental relative pressures up to high salt concentrations, despite the inability of these models to predict cohesive energies of the solutions or the salts. A thermodynamic analysis reveals that the coarse-grained models achieve the experimental activity coefficients of water in solution through a compensation of severely underestimated hydration and vaporization free energies of the salts. Our results suggest that coarse-grained models developed to replicate the hydration structure and the effective ion–ion attraction in solution may lead to this compensation. Moreover, they suggest an avenue for the design of coarse-grained models that accurately reproduce the activity coefficients of solutions.

Crystallization of ice from deeply supercooled water and amorphous ices – a process of fundamental importance in the atmosphere, interstellar space, and cryobiology – results in stacking disordered ices with a wide range of metastabilities with respect to hexagonal ice. The structural origin of this high variability, however, has not yet been elucidated. Here we use molecular dynamics simulations with the mW water model to characterize the structure of ice freshly grown from supercooled water at temperatures from 210 to 270 K, the thermodynamics of stacking faults, line defects, and interfaces, and to elucidate the interplay between kinetics and thermodynamics in determining the structure of ice. In agreement with experiments, the ice grown in the simulations is stacking disordered with a random distribution of cubic and hexagonal layers, and a cubicity that decreases with growth temperature. The former implies that the cubicity of ice is determined by processes at the ice/liquid interface, without memory of the structure of buried ice layers. The latter indicates that the probability of building a cubic layer at the interface decreases upon approaching the melting point of ice, which we attribute to a more efficient structural equilibration of ice at the liquid interface as the driving force for growth wanes. The free energy cost for creating a pair of cubic layers in ice is 8.0 J mol−1 in experiments, and 9.7 ± 1.9 J mol−1 for the mW water model. This not only validates the simulations, but also indicates that dispersion in cubicity is not sufficient to explain the large energetic variability of stacking disordered ices. We compute the free energy cost of stacking disorder, line defects, and interfaces in ice and conclude that a characterization of the density of these defects is required to predict the degree of metastability and vapor pressure of atmospheric ices.

Water, silicon, silica, and other liquids that favor tetrahedral order display thermodynamic, dynamic, and structural anomalies in the pressure range in which they form tetrahedrally coordinated crystals. The tetrahedral order in these liquids is induced by anisotropic hydrogen bonding or covalent interactions, or, in ionic melts, by an appropriate size ratio of the ions. Simple isotropic two-length scale models have been extensively used to understand the origin of anomalies in complex liquids. However, single-component isotropic liquids characterized to date generally do not stabilize tetrahedral crystals, and in the few cases that they do, it was found that the liquids do not display anomalies in the region of the tetrahedral crystal. This poses the question of whether it is possible for isotropic pair potentials to display water-like phase behavior and anomalies. In this work, we use molecular dynamics simulations to investigate the phase behavior and the existence and loci of anomalies of a single-component purely repulsive isotropic pair potential that stabilizes diamond in the ground state over a wide range of pressures. We demonstrate that, akin to water, silica, and silicon, the isotropic potential of Marcotte, Stillinger, and Torquato (MST) presents structural, dynamic, and thermodynamic anomalies in the region of stability of the tetrahedral crystal. The regions of anomalies of MST are nested in the T–p plane following the same hierarchy as in silica: the region of diffusional anomalies encloses the region of structural anomalies, which in turn contains the region of thermodynamic anomalies. To our knowledge, MST is the first example of pair potential for which water-like anomalies are associated with the formation of tetrahedral order.

The morphology of liquid–liquid phase separated aerosols has a strong impact on their rate of gas and water uptake, and the type and rate of heterogeneous reactions in the atmosphere. However, it is extremely challenging to experimentally distinguish different morphologies (core–shell or partial wetting) of aerosols and to quantify the extent of wetting between the two phases. The aim of this work is to quantitatively predict the morphology of liquid–liquid aerosols from fundamental physical properties of the aerosol phases. We determine the equilibrium structure of liquid–liquid phase separated aerosols through free energy minimization and predict that the contact angle between the two liquids in the aerosol depends on the composition but not the amount of each phase. We demonstrate that for aerosols of diameter larger than ∼100 nm, the equilibrium contact angle can be accurately predicted from the surface tensions of each liquid with the vapor and between the two liquids through an expression that is identical to Young’s equation. The internal structure of smaller, ultrafine aerosols depends also on the value of the line tension between the two liquids and the vapor. The thermodynamic model accurately predicts the experimental morphology, core–shell or partial wetting, of all aerosols for which surface tensions are provided in the literature, and provides contact angles that cannot be accurately determined with state of the art experimental methods. We find that the contact angle of model atmospheric aerosols is rarely higher than 30°. We validate the thermodynamic predictions through molecular simulations of nonane–water droplets, and use the simulation data to compute line tension values that are in good agreement with theory and the analysis from experimental data in water–nonane droplets. Our finding of a simple analytical equation to compute the contact angle of liquid–liquid droplets should have broad application for the prediction of the morphology of two-phase atmospheric aerosols and its impact in their chemistry.

Clathrate hydrates are crystals in which water molecules form hydrogen-bonded cages that enclose small nonpolar molecules, such as methane. In the laboratory, clathrates are customarily synthesized from ice and gas guest under conditions for which homogeneous nucleation of hydrates is not possible. It is not known how ice assists in the nucleation of clathrate hydrates or how ice forms on clathrate hydrate in the case of self-preservation. There is no lattice matching between any plane of ice and clathrate hydrates; therefore, an interfacial transition layer has to form to connect the two crystals. Here, we use molecular dynamic simulations to study the structure of ice–clathrate interfaces produced by alignment and equilibration of the crystals, competitive growth of ice and clathrate from a common solution, nucleation of hydrate in the presence of a growing ice front, and nucleation of ice in the presence of clathrate hydrates. We find that the interfacial transition layer between ice and clathrate is always disordered and has a typical width of two to three water layers. Water in the interfacial transition layer has tetrahedral order lower than either ice or clathrate and higher than liquid water under the same thermodynamic conditions. The potential energy of the water in the interfacial transition layer is between those in liquid water and the crystals. Our results suggest that the disordered interfacial transition layer could assist in the heterogeneous nucleation of clathrates from ice and ice from clathrates by providing a lower surface free energy than the ice–liquid and clathrate–liquid interfaces.

Clathrate hydrates and ice I are the most abundant crystals of water. The study of their nucleation, growth, and decomposition using molecular simulations requires an accurate and efficient algorithm that distinguishes water molecules that belong to each of these crystals and the liquid phase. Existing algorithms identify ice or clathrates, but not both. This poses a challenge for cases in which ice and hydrate coexist, such as in the synthesis of clathrates from ice and the formation of ice from clathrates during self-preservation of methane hydrates. Here we present an efficient algorithm for the identification of clathrate hydrates, hexagonal ice, cubic ice, and liquid water in molecular simulations. CHILL+ uses the number of staggered and eclipsed water–water bonds to identify water molecules in cubic ice, hexagonal ice, and clathrate hydrate. CHILL+ is an extension of CHILL (Moore et al. Phys. Chem. Chem. Phys. 2010, 12, 4124–4134), which identifies hexagonal and cubic ice but not clathrates. In addition to the identification of hydrates, CHILL+ significantly improves the detection of hexagonal ice up to its melting point. We validate the use of CHILL+ for the identification of stacking faults in ice and the nucleation and growth of clathrate hydrates. To our knowledge, this is the first algorithm that allows for the simultaneous identification of ice and clathrate hydrates, and it does so in a way that is competitive with respect to existing methods used to identify any of these crystals.

Atmospheric aerosols have a strong influence on Earth’s climate. Elucidating the physical state and internal structure of atmospheric aqueous aerosols is essential to predict their gas and water uptake, and the locus and rate of atmospherically important heterogeneous reactions. Ultrafine aerosols with sizes between 3 and 15 nm have been detected in large numbers in the troposphere and tropopause. Nanoscopic aerosols arising from bubble bursting of natural and artificial seawater have been identified in laboratory and field experiments. The internal structure and phase state of these aerosols, however, cannot yet be determined in experiments. Here we use molecular simulations to investigate the phase behavior and internal structure of liquid, vitrified, and crystallized water–salt ultrafine aerosols with radii from 2.5 to 9.5 nm and with up to 10% moles of ions. We find that both ice crystallization and vitrification of the nanodroplets lead to demixing of pure water from the solutions. Vitrification of aqueous nanodroplets yields nanodomains of pure low-density amorphous ice in coexistence with vitrified solute rich aqueous glass. The melting temperature of ice in the aerosols decreases monotonically with an increase of solute fraction and decrease of radius. The simulations reveal that nucleation of ice occurs homogeneously at the subsurface of the water–salt nanoparticles. Subsequent ice growth yields phase-segregated, internally mixed, aerosols with two phases in equilibrium: a concentrated water–salt amorphous mixture and a spherical cap-like ice nanophase. The surface of the crystallized aerosols is heterogeneous, with ice and solution exposed to the vapor. Free energy calculations indicate that as the concentration of salt in the particles, the advance of the crystallization, or the size of the particles increase, the stability of the spherical cap structure increases with respect to the alternative structure in which a core of ice is fully surrounded by solution. We predict that micrometer-sized particles and nanoparticles have the same equilibrium internal structure. The variation of liquid–vapor surface tension with solute concentration is a key factor in determining whether a solution-embedded ice core or vapor-exposed ice cap is the equilibrium structure of the aerosols. In agreement with experiments, we predict that the structure of mixed-phase HNO3–water particles, representative of polar stratospheric clouds, consists of an ice core surrounded by freeze-concentrated solution. The results of this work are important to determine the phase state and internal structure of sea spray ultrafine aerosols and other mixed-phase particles under atmospherically relevant conditions.

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 mXREM 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 mXREM 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 mXREM models and liquid–ice equilibrium properties show reasonable agreement with the results of the corresponding atomistic models. The mXREM 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/2005REM 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.

Carbonaceous particles account for 10% of the particulate matter in the atmosphere. Atmospheric oxidation and aging of soot modulates its ice nucleation ability. It has been suggested that an increase in the ice nucleation ability of aged soot results from an increase in the hydrophilicity of the surfaces upon oxidation. Oxidation, however, also impacts the nanostructure of soot, making it difficult to assess the separate effects of soot nanostructure and hydrophilicity in experiments. Here we use molecular dynamics simulations to investigate the effect of changes in hydrophilicity of model graphitic surfaces on the freezing temperature of ice. Our results indicate that the hydrophilicity of the surface is not in general a good predictor of ice nucleation ability. We find a correlation between the ability of a surface to promote nucleation of ice and the layering of liquid water at the surface. The results of this work suggest that ordering of liquid water in contact with the surface plays an important role in the heterogeneous ice nucleation mechanism.

Water flow coupled to the migration of ions is an important aspect of the performance of polymer electrolyte membrane fuel cells. The water gradients arising from the operation of fuel cells can result in flooding and drying-out of the electrodes and drying of regions of the membrane, with concomitant losses in conductivity and efficiency. The electro-osmotic drag coefficient measures the ratio between the flow of solvent molecules to that of a charged species toward an electrode in the presence of an applied electric field. The effects of variables such as pore radius, crowding, temperature, electric field strength, and ion concentration on the mobility of ions and accompanying water molecules in an applied electric field are still not well understood. Here, we investigate these factors with coarse-grained molecular simulations using an efficient model of water and sodium chloride ions and compare these results with those from previous experiments on proton exchange membranes as well as new experimental results for an anion exchange membrane. The anion exchange membranes have a smaller value of Kdrag than the proton-exchange membranes, which may be attributed to smaller water domains and a different charge carrier (hydroxide instead of protons). We directly determine the role of pore size on Kdrag and confirm that narrower pores result in less electro-osmotic drag. Our simulations show that Kdrag is sensitive to the interaction of the charge carrier with water molecules. The results of this work suggest that the most promising approach to minimize electro-osmotic drag while maintaining adequate ion conductivity is to control the morphology of the membrane structure at the microscopic level.

Coarse-grained models are becoming increasingly popular due to their ability to access time and length scales that are prohibitively expensive with atomistic models. However, as a result of decreasing the degrees of freedom, coarse-grained models often have diminished accuracy, representability, and transferability compared with their finer grained counterparts. Uncertainty quantification (UQ) can help alleviate this challenge by providing an efficient and accurate method to evaluate the effect of model parameters on the properties of the system. This method is useful in finding parameter sets that fit the model to several experimental properties simultaneously. In this work we use UQ as a tool for the evaluation and optimization of a coarse-grained model. We efficiently sample the five-dimensional parameter space of the coarse-grained monatomic water (mW) model to determine what parameter sets best reproduce experimental thermodynamic, structural and dynamical properties of water. Generalized polynomial chaos (gPC) was used to reconstruct the analytical surfaces of density, enthalpy of vaporization, radial and angular distribution functions, and diffusivity of liquid water as a function of the input parameters. With these surfaces, we evaluated the sensitivity of these properties to perturbations of the model input parameters and the accuracy and representability of the coarse-grained models. In particular, we investigated what is the optimum length scale of the water–water interactions needed to reproduce the properties of liquid water with a monatomic model with two- and three-body interactions. We found that there is an optimum cutoff length of 4.3 Å, barely longer than the size of the first neighbor shell in water. As cutoffs deviate from this optimum value, the ability of the model to simultaneously reproduce the structure and thermodynamics is severely diminished.

Clathrate hydrates are nonstoichiometric compounds comprised of a hydrogen-bonded water network that forms polyhedral cages that can be occupied by small guest molecules. Clathrates are candidate materials for storage and transportation of methane and H2. Promoter molecules, such as THF, reduce the pressure or temperature needed to form clathrates of these gases, resulting in the formation of binary clathrates with the promoter molecule hosted in the large cages of the hydrate. In this work, we study the growth and occupancy of binary clathrates as a function of supercooling of the solution using molecular dynamics simulations with the mW water model and small and large guest molecules of sizes similar to those of H2 and THF, respectively, but that are both highly soluble in water and produce single hydrates with identical melting point. The large guest molecules only fit into the 51264 cages, while the small guest molecules can fit into both types of cages. We find that the large guest act as a kinetic promoter for growth, increasing the rate of uptake of small guests into the clathrate. Our results also indicate that the growth of binary clathrates is limited by the arrangement of guest molecules in the large 51264 cages at the clathrate/solution interface. The occupancy of large cages of binary clathrates can be tuned by varying the growth temperature. The simulations indicate that with increasing supercooling there is an increase in the percentage of 51264 cages occupied by the small guest molecules at the expense of the large guest molecules, while the occupancy of 512 cages remains relatively constant. The results of this work show that the composition of clathrates grown at high driving force does not necessarily reflect the composition of the most stable phase.

Clathrate hydrates are crystals in which water forms a network of fully hydrogen-bonded polyhedral cages that contain small guests. Clathrate hydrates occur mostly in two cubic crystal polymorphs, sI and sII. Bromine is one of two guests that yield a hydrate with the tetragonal structure (TS), the topological dual of the Frank–Kasper σ phase. There has been a long-standing disagreement on whether bromine hydrate also forms metastable sI and sII crystals. To date there are no data on the thermodynamic range of stability (e.g., the melting temperatures) of the metastable polymorphs. Here we use molecular dynamics simulations with the coarse-grained model of water mW to (i) investigate the thermodynamic stability of the empty and guest-filled the sI, sII, TS, and HS-I hydrate polymorphs, (ii) develop a coarse-grained model of bromine compatible with mW water, and (iii) evaluate the stability of the bromine hydrate polymorphs. The mW model predicts the same relative energy of the empty clathrate polymorphs and the same phase diagram as a function of water–guest interaction than the fully atomistic TIP4P water model. There is a narrow region in water-guest parameter space for which TS is marginally more stable than sI or sII. We parametrize a coarse-grained model of bromine compatible with mW water and use it to determine the order of stability of the bromine hydrate polymorphs. The melting temperatures of the bromine hydrate polymorphs predicted by the coarse-grained model are 281 ± 1 K for TS, 279 ± 1 K for sII, and 276 ± 1 K for sI. The closeness of the melting temperatures supports the plausibility of formation of metastable sII and sI bromine hydrates.

Water nanoparticles play an important role in atmospheric processes, yet their equilibrium and nonequilibrium liquid–ice phase transitions and the structures they form on freezing are not yet fully elucidated. Here we use molecular dynamics simulations with the mW water model to investigate the nonequilibrium freezing and equilibrium melting of water nanoparticles with radii R between 1 and 4.7 nm and the structure of the ice formed by crystallization at temperatures between 150 and 200 K. The ice crystallized in the particles is a hybrid form of ice I with stacked layers of the cubic and hexagonal ice polymorphs in a ratio approximately 2:1. The ratio of cubic ice to hexagonal ice is insensitive to the radius of the water particle and is comparable to that found in simulations of bulk water around the same temperature. Heating frozen particles that contain multiple crystallites leads to Ostwald ripening and annealing of the ice structures, accompanied by an increase in the amount of ice at the expense of the liquid water, before the particles finally melt from the hybrid ice I to liquid, without a transition to hexagonal ice. The melting temperatures Tm of the nanoparticles are not affected by the ratio of cubic to hexagonal layers in the crystal. Tm of the ice particles decreases from 255 to 170 K with the particle size and is well described by the Gibbs–Thomson equation, Tm(R) = Tmbulk – KGT/(R – d), with constant KGT = 82 ± 5 K·nm and a premelted liquid of width d = 0.26 ± 0.05 nm, about one monolayer. The freezing temperatures also decrease with the particles’ radii. These results are important for understanding the composition, freezing, and melting properties of ice and liquid water particles under atmospheric conditions.

The role of water (thermo)dynamics is crucial in molecular recognition and self-assembly. Here, we study a prototype cavity–ligand system as a model for hydrophobic noncovalent binding. Two alternative molecular dynamics simulation resolutions are employed and the resulting structural, dynamic, and thermodynamic properties compared: first, a coarse-grained (CG) resolution based on the previously reported and validated methane-like M solute and mW water models; second, an atomic-level (AL) resolution based on the popular OPLS united atom methane and the TIP4P water models. The CG model reproduces, as a function of the cavity–ligand distance, (1) the water occupancy of the cavity, (2) the cavity–ligand potential of mean force (free energy) and its temperature dependence, and (3) some of the major qualitative features of the thermodynamic signatures (free energy, enthalpy, and entropy) for cavity–ligand association of the AL model. The limits of the CG and AL models in this context are also discussed with comparison to experimental data. Our study suggests that CG simulation with models that include the translational contribution of water and anisotropic “hydrogen-bond”-like interactions could reproduce the thermodynamics of molecular recognition and water-driven assembly in complex macromolecular systems and nanoscale processes with convenient computational time savings.

A premelted layer of water wets the surface of ice at temperatures below the melting temperature. Experiments suggest that this quasi-liquid layer may play an important role in the nucleation of clathrate hydrates from ice. Nevertheless, the structure of the quasi-liquid layer of ice in the presence of methane or other clathrate-forming gases has not yet been elucidated. In this work, we perform large-scale molecular dynamic simulations with a coarse-grained molecular model to investigate the properties of the quasi-liquid layer of ice in the presence of methane gas under pressure. We characterize the structure and thickness of the ice/methane and ice/vacuum interfaces, and the solubility of methane in the premelted layer as a function of temperature. We find that the width of the quasi-liquid layer fluctuates between 5 and 45 Å in the presence of a methane-like solute at temperatures within 1 K of the melting point. The width of the quasi-liquid layer of ice at temperatures lower than 270 K is less than the diameter of a water dodecahedron, the smallest cage that constitutes the clathrates. The simulations indicate that, when the premelting layer is wider than 10 Å, the structure of water and solubility of methane in the center of the quasi-liquid layer are the same as in bulk liquid water at the same temperature. These results are relevant for understanding the mechanism of formation of methane hydrate clathrates from ice.

The liquid–vapor transition in cylindrical pores is studied as a function of pore size and hydrophilicity through molecular dynamics simulations with the mW coarse-grained model of water. We identify two distinct filling mechanisms, depending on whether the water–pore interaction is smaller or larger than the water–water interaction. In the former case (that we term hydrophobic pore), the formation of the condensed phase proceeds gradually with filling, through the nucleation of a water cluster which grows toward the center of the cavity. In hydrophilic pores, instead, the condensed phase develops in conditions of supersaturation, which in principle become more extreme with increasing pore radius and surface affinity. For highly hydrophilic interfaces (those with adsorption energy for water above 10 kcal/mol), the equilibrium and dynamical properties of water in confinement turn out to be practically independent of water affinity.

We use molecular dynamics simulations to investigate the coexistence between confined ice and liquid water as a function of temperature for a series of cylindrical nanopores with water–wall interactions ranging from strongly hydrophilic to very hydrophobic. In agreement with previous results from experiments, we find that the ice formed in the nanopores is a hybrid ice I with stacks of cubic and hexagonal layers and that the melting temperature of the nanoconfined ice is strongly dependent on the radius of the pore but rather insensitive to the hydrophilicity of the pore surface. We find a premelted liquid layer in coexistence with the confined ice down to the lowest temperatures of this study, 50 K below the melting temperatures of the confined ices. The fraction of water in the premelted liquid layer decreases with increasing hydrophobicity of the pore wall, but it does not vanish even for the most hydrophobic nanopores. The simulations suggest that the decrease in the fraction of water in the liquid layer with increasing hydrophobicity corresponds partly to a decrease in its width but also to a decrease in its effective density: the premelted liquid layer becomes sparser on decreasing the water–wall attraction. Our results indicate that agreement in the melting temperatures of water nanopores functionalized with different moieties does not imply identical fraction of nonfreezable water in these materials.

Clathrate hydrates mostly occur in two cubic crystal structures, sI and sII. Cross-nucleation between these clathrate crystals has been observed in simulations and may be relevant to the transformation between clathrate polymorphs reported in experiments. Nevertheless, the mechanism by which clathrate crystals cross-nucleate and the structure of the interface between the distinct crystals have not yet been fully characterized. In this work, we use extensive molecular dynamics simulations to investigate the structure of the clathrate/solution interface for sI and sII guest-free and methane-filled hydrates at different degrees of supercooling and the mechanism of cross-nucleation between clathrate polymorphs. We find that 51263 water cages, which are not native to the sI or sII crystals, occur assiduously in their interfaces with the solution and play a central role in the mechanism of cross-nucleation of clathrate hydrates: cross-nucleation between sI and sII requires the formation of an interfacial layer tiled by 51263 cages connected by dodecahedra. We characterize the structure of the interfacial layer, estimate the size of the critical surface nucleus required for its formation, and assess the role of other variables in the reaction coordinate of cross-nucleation of clathrate hydrates. In agreement with previous reports of cross-nucleation of quite different systems, we observe cross-nucleation of clathrate hydrates both from the stable to the metastable crystal and from the metastable to the stable hydrate. The new crystal that forms is, in all cases, the one that has the fastest growth rate.

Recent studies reveal that amorphous intermediates are involved in the formation of clathrate hydrates under conditions of high driving force, raising two questions: first, how could amorphous nuclei grow into crystalline clathrates and, second, whether amorphous nuclei are intermediates in the formation of clathrate crystals for temperatures close to equilibrium. In this work, we address these two questions through large-scale molecular simulations. We investigate the stability and growth of amorphous and crystalline clathrate nuclei and assess the thermodynamics and kinetic factors that affect the crystallization pathway of clathrates. Our calculations show that the dissociation temperature of amorphous clathrates is just 10% lower than for the crystals, facilitating the formation of metastable amorphous intermediates. We find that, at any temperatures, the critical crystalline nuclei are smaller than critical amorphous nuclei. The temperature dependence of the critical nucleus size is well described by the Gibbs−Thomson relation, from which we extract a liquid-crystal surface tension in excellent agreement with experiments. Our analysis suggests that at high driving force the amorphous nuclei may be kinetically favored over crystalline nuclei because of lower free energy barriers of formation. We investigated the role of the initial structure and size of the nucleus on the subsequent growth of the clathrates, and found that both amorphous and sI crystalline nuclei yield crystalline clathrates. Interestingly, growth of the metastable sII crystal polymorph is always favored over the most stable sI crystal, revealing kinetic control of the growth and indicating that a further step of ripening from sII to sI is needed to reach the stable crystal phase. The latter results are in agreement with the observed metastable formation of sII CO2 and CH4 clathrate hydrates and their slow conversion to sI under experimental conditions.

Ice crystallized below 200 K has the diffraction pattern of a faulty cubic ice, and not of the most stable hexagonal ice polymorph. The origin and structure of this faulty cubic ice, presumed to form in the atmosphere, has long been a puzzle. Here we use large-scale molecular dynamics simulations with the mW water model to investigate the crystallization of water at 180 K and elucidate the development of cubic and hexagonal features in ice as it nucleates, grows and consolidates into crystallites with characteristic dimensions of a few nanometres. The simulations indicate that the ice crystallized at 180 K contains layers of cubic ice and hexagonal ice in a ratio of approximately 2 to 1. The stacks of hexagonal ice are very short, mostly one and two layers, and their frequency does not seem to follow a regular pattern. In spite of the high fraction of hexagonal layers, the diffraction pattern of the crystals is, as in the experiments, almost identical to that of cubic ice. Stacking of cubic and hexagonal layers is observed for ice nuclei with as little as 200 water molecules, but a preference for cubic ice is already well developed in ice nuclei one order of magnitude smaller: the critical ice nuclei at 180 K contain approximately ten water molecules in their core and are already rich in cubic ice. The energies of the cubic-rich and hexagonal-rich nuclei are indistinguishable, suggesting that the enrichment in cubic ice does not have a thermodynamic origin.

Water confined in narrow nanopores that prevent ice crystallization is usually studied as a means to understand the anomalous behavior of bulk liquid water. Nevertheless, there is no agreement on the similarity of the thermodynamics of bulk and nanoconfined liquid water. In this work, we use molecular dynamics simulations with the mW water model to investigate the phase behavior of liquid water in bulk and confined in a 1.5 nm cylindrical pore with water–surface interactions identical to water–water interactions. Through analysis of the isochors of bulk liquid water V(T,p) we extrapolate the locus of a putative liquid–liquid critical point (LLCP) in bulk mW water at 190 K and 1215 atm. This is a “virtual LLCP”, as it would lie in a region of the phase diagram where fast crystallization of water impedes equilibration of the liquid. We find that confinement has a weak effect on the loci of the thermodynamic anomalies: the maxima in density and heat capacity of confined water occur at T,p similar to that in bulk. The heat capacity peak of confined water is due to a transformation within the confined liquid; we verify that ice does not form in the pores. Confined mW water presents a heat capacity maximum Cpmax(p) up to the highest pressure we investigated, 4000 atm. The magnitude of the heat capacity peak Cpmax(p) has a nonmonotonous dependence with pressure, attaining a maximum at conditions close to those of the locus of the bulk water’s virtual LLCP. We do not find, however, direct evidence of a first-order liquid–liquid transition in confined water for pressures above or below the locus of the maximum response function. The extreme value of the response functions of confined water could be a rounded manifestation of an equivalent feature in the free energy surface of bulk water.

We investigate the melting and formation of ice in partially filled hydrophilic and hydrophobic nanopores of 3 nm diameter using molecular dynamics simulations with the mW water model. Above the melting temperature, the partially filled nanopores contain two water phases in coexistence: a condensed liquid plug and a surface-adsorbed phase. It has been long debated in the literature whether the surface-adsorbed phase is involved in the crystallization. We find that only the liquid plug crystallizes on cooling, producing ice I with stacks of hexagonal and cubic layers. The confined ice is wetted by a premelted liquid layer that persists in equilibrium with ice down to temperatures well below its melting point. The liquid–ice transition is first-order-like but rounded. We determine the temperature and enthalpy of melting as a function of the filling fraction of the pore. In agreement with experiments, we find that the melting temperature of the nanoconfined ice is strongly depressed with respect to the bulk Tm, it depends weakly on the filling fraction and is insensitive to the hydrophobicity of the pore wall. The state of water in the crystallized hydrophilic and hydrophobic pores, however, is not the same: the hydrophobic pore has a negligible density of the surface-adsorbed phase and higher fraction of water in the ice phase than the hydrophilic pore. The widths of the ice cores are nevertheless comparable for the hydrophobic and hydrophilic pores, and this may explain their almost identical melting temperatures. The enthalpy of melting ΔHm, when normalized by the actual amount of ice in the pore, is indistinguishable for the hydrophobic and hydrophilic pores, insensitive to the filling fraction, and within the error bars, the same as the difference in enthalpy between bulk liquid and bulk ice evaluated at the temperature of melting of ice in the nanopores.

We use large-scale molecular dynamics simulations to investigate the phase transformation of aqueous solutions of electrolytes cooled at the critical rate to avoid the crystallization of ice. Homogeneous liquid solutions with up to 20% moles of ions demix on cooling producing nanophase segregated glasses with characteristic dimensions of phase segregation of about 5 nm. The immiscibility is driven by the transformation of water to form a four-coordinated low-density liquid (LDL) as it crosses the liquid−liquid transformation temperature TLL of the solution. The ions cannot be incorporated into the tetrahedral LDL network and are expelled to form a solute-rich water nanophase. The simulations quantitatively reproduce the relative amounts of low and high-density liquid water as a function of solute content in LiCl glasses [Suzuki and Mishima, Phys. Rev. Lett.2000, 85, 1322−1325] and provide direct evidence of segregation in aqueous glasses and their dimensions of phase segregation.

Solvation by water and ions has been shown to be vitally important for biological molecules, yet fully atomistic simulations of large biomolecules remain a challenge due to their high computational cost. The effect of solvation is the most pronounced in polyelectrolytes, of which DNA is a paradigmatic example. Coarse-grained (CG) representations have been developed to model the essential physics of the DNA molecule, yet almost without exception, these models replace the water and ions by implicit solvation in order to significantly reduce the computational expense. This work introduces the first coarse-grained model of DNA solvated explicitly with water and ions. To this end, we combined two established CG models; the recently developed mW-ion model [DeMille, R. C.; Molinero, V. J. Chem. Phys. 2009, 131, 034107], which reproduces the structure of aqueous ionic solutions without electrostatic interactions, was coupled to the three-sites-per-nucleotide (3SPN) CG model of DNA [Knotts, T. A., IV; et al. J. Chem. Phys. 2007, 126, 084901]. Using atomistic simulations of d(CGCGAATTCGCG)2 as a reference, we optimized the coarse-grained interactions between DNA and solvent to reproduce the solvation structure of water and ions around CG DNA. The resulting coarse-grained model of DNA explicitly solvated by ions and water (mW/3SPN-DNA) exhibits base-pair specificity and ion-condensation effects and it is 2 orders of magnitude computationally more efficient than atomistic models. We describe the parametrization strategy and offer insight into how other CG models may be combined with a coarse-grained solvent model such as mW-ion.

Quasicrystals are structures with long-range order and no translational periodicity. Monatomic quasicrystals were predicted for model potentials, but no single-component atomic quasicrystal of an actual element has been reported to date. A dodecagonal quasicrystal was recently predicted to form in bilayer water. Water and silicon present striking similarities in their phase behavior, raising the question of whether quasicrystals may occur in silicon. Here, we show, using molecular simulations, that a confined silicon bilayer forms a quasicrystal upon compression between smooth surfaces. The quasicrystal is stable in a narrow region of the phase diagram and forms spontaneously upon cooling the liquid bilayer in a wide range of pressures. Cooling the liquid between atomically detailed plates incommensurate with the quasicrystal leads to its spontaneous formation at 1 atm of lateral pressure. This suggests that the silicon quasicrystal could be obtained in experiments at room pressure by tuning the structure and interactions of the surfaces.Keywords: confined; crystallization; interfacial; phase transitions; quasicrystal; silicon; water;

The nucleation, growth, structure and melting of ice in 3 nm diameter hydrophilic nanopores are studied through molecular dynamics simulations with the mW water model. The melting temperature of water in the pore was Tporem = 223 K, 51 K lower than the melting point of bulk water in the model and in excellent agreement with experimental determinations for 3 nm silica pores. Liquid and ice coexist in equilibrium at the melting point and down to temperatures as low as 180 K. Liquid water is located at the interface of the pore wall, increasing from one monolayer at the freezing temperature, Tporef = 195 K, to two monolayers a few degrees below Tporem. Crystallization of ice in the pore occurs through homogeneous nucleation. At the freezing temperature, the critical nucleus contains ∼75 to 100 molecules, with a radius of gyration similar to the radius of the pore. The critical nuclei contain features of both cubic and hexagonal ice, although stacking of hexagonal and cubic layers is not defined until the nuclei reach ∼150 molecules. The structure of the confined ice is rich in stacking faults, in agreement with the interpretation of X-ray and neutron diffraction experiments. Though the presence of cubic layers is twice as prevalent as hexagonal ones, the crystals should not be considered defective Ic as sequences with more than three adjacent cubic (or hexagonal) layers are extremely rare in the confined ice.

Methane is the prototypic hydrophobic molecule; it has an extremely low solubility in liquid water that leads to phase segregation. On the other hand, at moderate pressures and room temperature, water and methane form hydrate clathrate crystals with a methane to water ratio up to a 1000 times higher than the saturated aqueous phase. This apparent dichotomy points to a subtle balance between the strong water−water hydrogen bonding, responsible for the hydrophobic effect, and water−methane attraction. Capturing these nuances with molecular models requires an appropriate balance of intermolecular interactions. Here we present such a coarse-grained molecular model of water and methane that represents each molecule by a single particle interacting through very short-range interaction potentials. The model is based on the monatomic model of water mW [Molinero, V.; Moore, E. B. J. Phys. Chem. B 2009, 113, 4008] and is between 2 and 3 orders of magnitude more computationally efficient than atomistic models with Ewald sums. The coarse-grained model of this study reproduces the solubility and hydration number of methane in liquid water, the surface tension of the water−methane interface and the equilibrium melting temperature of methane hydrate clathrates with structures sI and sII. To the best of our knowledge this is the first force-field, atomistic or coarse-grained, that reproduces these range of properties of liquid and solid phases of water and methane, making it an efficient and accurate model for the study of the mechanisms of nucleation and growth of clathrates. We expect that the results of this work will also be useful for the modeling of the hydrophobic assembly in aqueous solutions and the development of coarse-grained models of biomolecules with explicit solvation.

Understanding the microscopic mechanism of nucleation of clathrate hydrates is important for their use in hydrogen storage, CO2 sequestration, storage and transport of natural gas, and the prevention of the formation of hydrate plugs in oil and gas pipelines. These applications involve hydrate guests of varied sizes and solubility in water that form different hydrate crystal structures. Nevertheless, molecular studies of the mechanism of nucleation of hydrates have focused on the single class of small hydrophobic guests that stabilize the sI crystal. In this work, we use molecular dynamics simulations with a very efficient coarse-grained model to elucidate the mechanisms of nucleation of clathrate hydrates of four model guests that span a 2 orders of magnitude range in solubility in water and that encompass sizes which stabilize each one a different hydrate structure (sI and sII, with and without occupancy of the dodecahedral cages). We find that the overall mechanism of clathrate nucleation is similar for all guests and involves a first step of formation of blobs, dense clusters of solvent-separated guest molecules that are the birthplace of the clathrate cages. Blobs of hydrophobic guests are rarer and longer-lived than those for soluble guests. For each guest, we find multiple competing channels to form the critical nuclei, filled dodecahedral (512) cages, empty 512 cages, and a variety of filled large (5126n with n = 2, 3, and 4) clathrate cages. Formation of empty dodecahedra is an important nucleation channel for all but the smallest guest. The empty 512 cages are stabilized by the presence of guests from the blob in their first solvation shell. Under conditions of high supercooling, the structure of the critical and subcritical nuclei is mainly determined by the size of the guest and does not reflect the cage composition or ordering of the stable or metastable clathrate crystals.

We use extensive molecular dynamics simulations with the monatomic model of water (mW) to characterize the thermodynamics and kinetics of the liquid−vapor (wetting−drying) equilibrium of water confined between nanoscopic hydrophobic plates. The transition in confined water is first-order-like, with two well-defined states (wet and dry) separated by a free energy barrier. Different from its bulk counterpart, the confined system oscillates between liquid and vapor: the two phases coexist in time but not in space. Also different from the phase behavior in bulk, there is a finite range of the thermodynamic variables (e.g., temperature or separation between the plates) for which the liquid and vapor state coexist in dynamical equilibrium. We determine the range of temperatures and plate separations for which reversible oscillations can be observed between a stable and metastable phase, compute the time scales of the phase transition along the equilibrium coexistence line, and investigate the pathway for drying along simple collective coordinates that describe the opening of a vapor bubble. The results of the simulations are compared with a simple capillary model for the thermodynamics and transition state theory for the kinetics of phase oscillations.

Water and silicon are chemically dissimilar substances with common physical properties. Their liquids display a temperature of maximum density, increased diffusivity on compression, and they form tetrahedral crystals and tetrahedral amorphous phases. The common feature to water, silicon, and carbon is the formation of tetrahedrally coordinated units. We exploit these similarities to develop a coarse-grained model of water (mW) that is essentially an atom with tetrahedrality intermediate between carbon and silicon. mW mimics the hydrogen-bonded structure of water through the introduction of a nonbond angular dependent term that encourages tetrahedral configurations. The model departs from the prevailing paradigm in water modeling: the use of long-ranged forces (electrostatics) to produce short-ranged (hydrogen-bonded) structure. mW has only short-range interactions yet it reproduces the energetics, density and structure of liquid water, and its anomalies and phase transitions with comparable or better accuracy than the most popular atomistic models of water, at less than 1% of the computational cost. We conclude that it is not the nature of the interactions but the connectivity of the molecules that determines the structural and thermodynamic behavior of water. The speedup in computing time provided by mW makes it particularly useful for the study of slow processes in deeply supercooled water, the mechanism of ice nucleation, wetting-drying transitions, and as a realistic water model for coarse-grained simulations of biomolecules and complex materials.

We use molecular dynamics simulations with the monatomic water (mW) model to investigate the phase diagram, metastability, and growth of guest-free water clathrates of structure sI and sII. At 1 atm pressure, the simulated guest-free water clathrates are metastable with respect to ice and stable with respect to the liquid up to their melting temperatures, 245 ± 2 and 252 ± 2 K for sI and sII, respectively. We characterize the growth of the sI and sII water crystals from supercooled water and find that the clathrates are unable to nucleate ice, the stable crystal. We computed the phase relations of ice, guest-free sII clathrate, and liquid water from −1500 to 500 atm. The resulting phase diagram indicates that empty sII may be the stable phase of water at pressures lower than approximately −1300 atm and temperatures lower than 275 K. The simulations show that, even in the region of stability of the empty clathrates, supercooled liquid water crystallizes to ice. This suggests that the barrier for nucleation of ice is lower than that for clathrates. We find no evidence for the existence of the characteristic polyhedral clathrate cages in supercooled extended water. Our results show that clathrates do not need the presence of a guest molecule to grow, but they need the guest to nucleate from liquid water. We conclude that nucleation of empty clathrates from supercooled liquid water would be extremely challenging if not impossible to attain in experiments. We propose two strategies to produce empty water clathrates in laboratory experiments at low positive pressures.

The nucleation and growth of clathrate hydrates of a hydrophobic guest comparable to methane or carbon dioxide are studied by molecular dynamics simulations of two-phase systems. The crystallization proceeds in two steps: First, the guest molecules concentrate in “blobs”, amorphous clusters involving multiple guest molecules in water-mediated configurations. These blobs are in dynamic equilibrium with the dilute solution and give birth to the clathrate cages that eventually transform it into an amorphous clathrate nucleus. In a second step, the amorphous clathrate transforms into crystalline clathrate. At low temperatures, the system can be arrested in the metastable amorphous clathrate phase for times sufficiently long for it to appear as an intermediate in the crystallization of clathrates. The “blob mechanism” unveiled in this work synthesizes elements of the labile cluster and local structuring hypotheses of clathrate nucleation and bears strong analogies to the two-step mechanisms of crystallization of proteins and colloids.

Atmospheric aerosols can promote the heterogeneous nucleation of ice, impacting the radiative properties of clouds and Earth’s climate. The experimental investigation of heterogeneous freezing of water droplets by carbonaceous particles reveals widespread ice freezing temperatures. It is not known which structural and chemical characteristics of soot account for the variability in ice nucleation efficiency. Here we use molecular dynamics simulations to investigate the nucleation of ice from liquid water in contact with graphitic surfaces. We find that atomically flat carbon surfaces promote heterogeneous nucleation of ice, while molecularly rough surfaces with the same hydrophobicity do not. Graphitic surfaces and other surfaces that promote ice nucleation induce layering in the interfacial water, suggesting that the order imposed by the surface on liquid water may play an important role in the heterogeneous nucleation mechanism. We investigate a large set of graphitic surfaces of various dimensions and radii of curvature and find that variations in nanostructures alone could account for the spread in the freezing temperatures of ice on soot in experiments. We conclude that a characterization of the nanostructure of soot is needed to predict its ice nucleation efficiency.

The nucleation, growth, structure and melting of ice in 3 nm diameter hydrophilic nanopores are studied through molecular dynamics simulations with the mW water model. The melting temperature of water in the pore was Tporem = 223 K, 51 K lower than the melting point of bulk water in the model and in excellent agreement with experimental determinations for 3 nm silica pores. Liquid and ice coexist in equilibrium at the melting point and down to temperatures as low as 180 K. Liquid water is located at the interface of the pore wall, increasing from one monolayer at the freezing temperature, Tporef = 195 K, to two monolayers a few degrees below Tporem. Crystallization of ice in the pore occurs through homogeneous nucleation. At the freezing temperature, the critical nucleus contains ∼75 to 100 molecules, with a radius of gyration similar to the radius of the pore. The critical nuclei contain features of both cubic and hexagonal ice, although stacking of hexagonal and cubic layers is not defined until the nuclei reach ∼150 molecules. The structure of the confined ice is rich in stacking faults, in agreement with the interpretation of X-ray and neutron diffraction experiments. Though the presence of cubic layers is twice as prevalent as hexagonal ones, the crystals should not be considered defective Ic as sequences with more than three adjacent cubic (or hexagonal) layers are extremely rare in the confined ice.

The design of polymer electrolyte membranes with controlled water uptake is of high importance for high-performance fuel cells, because the water content of the membranes modulates their conductivity, chemical stability and mechanical strength. The water activity aw controls the equilibrium water uptake of a system. Predicting aw of materials is currently a daunting challenge for molecular simulations, because calculations of water activity require grand canonical simulations that are extremely expensive even with classical non-polarizable force fields. Moreover, force fields do not generally reproduce aw of solutions. Here, we first present a general strategy to parameterize force fields that reproduce the experimental aw of solutions, and then implement that strategy to re-parameterize the interactions in FFcomp, a coarse-grained model for hydrated polyphenylene oxide/trimethylamine chloride (PPO/TMACl) membranes in which the TMA cation is attached to the PPO backbone and the Cl anion is in the mobile water nanophase. Coarse-grained models based on short-ranged potentials successfully model fuel cell membranes and other concentrated aqueous electrolyte solutions because electrostatic interactions are highly screened in these systems. The new force field, FFpvap, differs from the original FFcomp only in the parameters of the ion–ion interactions, yet it reproduces aw in TMACl solutions with accuracy within 0.5 and 3% of the experimental value in all the concentration range relevant to the operation of fuel cell membranes. We find that the heat needed to vaporize water in solutions with as little as five water molecules per ion pair is essentially the same as in pure water, despite the strong water–ion interactions and their impact on the water activity. We review the literature to demonstrate that this is independent of the model and a general feature of water solutions. FFpvap reproduces the radial distribution functions and captures well the relative diffusivities of water and ions in the ionic solution as predicted by the reference atomistic GAFF-TIP4P/2005 model, while providing a hundred-fold gain in computing efficiency with respect to the atomistic model. With the backbone fragments inherited from FFcomp, the new FFpvap force field can be used to model hydrated polymer electrolyte membranes and advance the design of fuel cell membranes with controlled water uptake and conductivity.

Ice crystallized below 200 K has the diffraction pattern of a faulty cubic ice, and not of the most stable hexagonal ice polymorph. The origin and structure of this faulty cubic ice, presumed to form in the atmosphere, has long been a puzzle. Here we use large-scale molecular dynamics simulations with the mW water model to investigate the crystallization of water at 180 K and elucidate the development of cubic and hexagonal features in ice as it nucleates, grows and consolidates into crystallites with characteristic dimensions of a few nanometres. The simulations indicate that the ice crystallized at 180 K contains layers of cubic ice and hexagonal ice in a ratio of approximately 2 to 1. The stacks of hexagonal ice are very short, mostly one and two layers, and their frequency does not seem to follow a regular pattern. In spite of the high fraction of hexagonal layers, the diffraction pattern of the crystals is, as in the experiments, almost identical to that of cubic ice. Stacking of cubic and hexagonal layers is observed for ice nuclei with as little as 200 water molecules, but a preference for cubic ice is already well developed in ice nuclei one order of magnitude smaller: the critical ice nuclei at 180 K contain approximately ten water molecules in their core and are already rich in cubic ice. The energies of the cubic-rich and hexagonal-rich nuclei are indistinguishable, suggesting that the enrichment in cubic ice does not have a thermodynamic origin.

Crystallization of ice from deeply supercooled water and amorphous ices – a process of fundamental importance in the atmosphere, interstellar space, and cryobiology – results in stacking disordered ices with a wide range of metastabilities with respect to hexagonal ice. The structural origin of this high variability, however, has not yet been elucidated. Here we use molecular dynamics simulations with the mW water model to characterize the structure of ice freshly grown from supercooled water at temperatures from 210 to 270 K, the thermodynamics of stacking faults, line defects, and interfaces, and to elucidate the interplay between kinetics and thermodynamics in determining the structure of ice. In agreement with experiments, the ice grown in the simulations is stacking disordered with a random distribution of cubic and hexagonal layers, and a cubicity that decreases with growth temperature. The former implies that the cubicity of ice is determined by processes at the ice/liquid interface, without memory of the structure of buried ice layers. The latter indicates that the probability of building a cubic layer at the interface decreases upon approaching the melting point of ice, which we attribute to a more efficient structural equilibration of ice at the liquid interface as the driving force for growth wanes. The free energy cost for creating a pair of cubic layers in ice is 8.0 J mol−1 in experiments, and 9.7 ± 1.9 J mol−1 for the mW water model. This not only validates the simulations, but also indicates that dispersion in cubicity is not sufficient to explain the large energetic variability of stacking disordered ices. We compute the free energy cost of stacking disorder, line defects, and interfaces in ice and conclude that a characterization of the density of these defects is required to predict the degree of metastability and vapor pressure of atmospheric ices.