Water Freezing and Ice Melting

Comparing brute force to transition path sampling for gas hydrate nucleation with a flat interface: comments on time reversal symmetry

Fluid to solid nucleation is often investigated with the rare event method transition path sampling (TPS). I claim that the inherent irreversibility of solid nucleation, even at stationary conditions, calls into question TPS's applicability for determining solid nucleation mechanisms, especially for pre-critical behavior. Even when applied to a phenomenon which displays time reversal asymmetry like solid nucleation, TPS is a good means of exploring phase space and giving trends in post-critical structure, and its ability to facilitate nucleation rate and free energy calculations remains outstanding. Forward-only splitting and ratcheting methods such as forward flux sampling are more attractive for understanding nucleation mechanisms as they do not require time reversal symmetry, but at low driving forces may suffer from the same limitations as brute force: they may never make it to the first ratchet. Here I briefly summarize the TPS method and gas hydrate nucleation simulation literature, focusing on topics within both to facilitate a comparison of brute force hydrate nucleation to transition path sampling of hydrate nucleation. Perhaps anecdotally, the brute force technique results in more crystalline trajectories despite having higher driving forces than TPS. I maintain this difference is because of the inherent irreversibility of hydrate nucleation, meaning its pre-critical behavior cannot accurately be determined by the melting trajectories that comprise approximately half of the configurations in TPS's path ensemble.

Melting domain size and recrystallization dynamics of ice revealed by time-resolved x-ray scattering

The phase transition between water and ice is ubiquitous and one of the most important phenomena in nature. Here, we performed time-resolved x-ray scattering experiments capturing the melting and recrystallization dynamics of ice. The ultrafast heating of ice I is induced by an IR laser pulse and probed with an intense x-ray pulse which provided us with direct structural information on different length scales. From the wide-angle x-ray scattering (WAXS) patterns, the molten fraction, as well as the corresponding temperature at each delay, were determined. The small-angle x-ray scattering (SAXS) patterns, together with the information extracted from the WAXS analysis, provided the time-dependent change of the size and the number of liquid domains. The results show partial melting (~13%) and superheating of ice occurring at around 20 ns. After 100 ns, the average size of the liquid domains grows from about 2.5 nm to 4.5 nm by the coalescence of approximately six adjacent domains. Subsequently, we capture the recrystallization of the liquid domains, which occurs on microsecond timescales due to the cooling by heat dissipation and results to a decrease of the average liquid domain size.

Tracking cubic ice at molecular resolution

Ice is present everywhere on Earth and has an essential role in several areas, such as cloud physics, climate change and cryopreservation. The role of ice is determined by its formation behaviour and associated structure. However, these are not fully understood1. In particular, there is a long-standing debate about whether water can freeze to form cubic ice-a currently undescribed phase in the phase space of ordinary hexagonal ice2-6. The mainstream view inferred from a collection of laboratory data attributes this divergence to the inability to discern cubic ice from stacking-disordered ice-a mixture of cubic and hexagonal sequences7-11. Using cryogenic transmission electron microscopy combined with low-dose imaging, we show here the preferential nucleation of cubic ice at low-temperature interfaces, resulting in two types of separate crystallization of cubic ice and hexagonal ice from water vapour deposition at 102 K. Moreover, we identify a series of cubic-ice defects, including two types of stacking disorder, revealing the structure evolution dynamics supported by molecular dynamics simulations. The realization of direct, real-space imaging of ice formation and its dynamic behaviour at the molecular level provides an opportunity for ice research at the molecular level using transmission electron microscopy, which may be extended to other hydrogen-bonding crystals.

Simulating the nematic-isotropic phase transition of liquid crystal model via generalized replica-exchange method

The nematic-isotropic (NI) phase transition of 4-cyano-4'-pentylbiphenyl was simulated using the generalized replica-exchange method (gREM) based on molecular dynamics simulations. The effective temperature is introduced in the gREM, allowing for the enhanced sampling of configurations in the unstable region, which is intrinsic to the first-order phase transition. The sampling performance was analyzed with different system sizes and compared with that of the temperature replica-exchange method (tREM). It was observed that gREM is capable of sampling configurations at sufficient replica-exchange acceptance ratios even around the NI transition temperature. A bimodal distribution of the order parameter at the transition region was found, which is in agreement with the mean-field theory. In contrast, tREM is ineffective around the transition temperature owing to the potential energy gap between the nematic and isotropic phases.

The microscopic mechanism of bulk melting of ice

We study the initial stages of homogeneous melting of a hexagonal ice crystal at coexistence and at moderate superheating. Our trajectory-based computer simulation approach provides a comprehensive picture of the events that lead to melting, from the initial accumulation of 5+7 defects, via the formation of L-D and interstitial-vacancy pairs, to the formation of a liquid nucleus. Of the different types of defects that we observe to be involved in melting, a particular kind of 5+7 type defect (type 5) plays a prominent role as it often forms prior to the formation of the initial liquid nucleus and close to the site where the nucleus forms. Hence, like other solids, ice homogeneously melts via the prior accumulation of defects.

Simulation of fluid/gel phase equilibrium in lipid vesicles

Simulation of single component dipalmitoylphosphatidylcholine (DPPC) coarse-grained DRY-MARTINI lipid vesicles of diameter 10 nm (1350 lipids), 20 nm (5100 lipids) and 40 nm (17 600 lipids) is performed using statistical temperature molecular dynamics (STMD), to study finite size effects upon the order-disorder gel/fluid transition. STMD obtains enhanced sampling using a generalized ensemble, obtaining a flat energy distribution between upper and lower cutoffs, with little computational cost over canonical molecular dynamics. A single STMD trajectory of moderate length is sufficient to sample 20+ transition events, without trapping in the gel phase, and obtain well averaged properties. Phase transitions are analyzed via the energy-dependence of the statistical temperature, T_S(U). The transition temperature decreases with decreasing diameter, in agreement with experiment, and the transition changes from first order to borderline first-second order. The size- and layer-dependence of the structure of both stable phases, and of the pathway of the phase transition, are determined. It is argued that the finite size effects are primarily caused by the disruption of the gel packing by curvature. Inhomogeneous states with faceted gel patches connected by unusual fluid seams are observed at high curvature, with visually different structure in the inner and outer layers due to the different curvatures. Thus a simple physical picture describes phase transitions in nanoscale finite systems far from the thermodynamic limit.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

The accretion of the new ice layer on the surface of hexagonal ice crystal and the influence of the local electric field on this process

The process of creation of a new layer of ice on the basal plane and on the prism plane of a hexagonal ice crystal is analyzed. It is demonstrated that the ordering of water molecules in the already existing crystal affects the freezing. On the basal plane, when the orientations of water molecules in the ice block are random, the arrangement of the new layer in a cubic manner is observed more frequently-approximately 1.7 times more often than in a hexagonal manner. When the water molecules in the ice block are more ordered, it results in the predominance of the oxygen atoms or the hydrogen atoms on the most outer part of the surface of the ice block. In this case, the hexagonal structure is formed more frequently when the supercooling of water exceeds 10 K. This phenomenon is explained by the influence of the oriented electric field, present as a consequence of the ordering of the dipoles of water molecules in the ice block. This field modifies the structure of solvation water (i.e., the layer of water in the immediate vicinity of the ice surface). We showed that the structure of solvation water predetermines the kind of the newly created layer of ice. This effect is temperature-dependent: when the temperature draws nearer to the melting point, the cubic structure becomes the prevailing form. The temperature at which the cubic and the hexagonal structures are formed with the same probabilities is equal to about 260 K. In the case of the prism plane, the new layer that is formed is always the hexagonal one, which is independent of the arrangement of water molecules in the ice block and is in agreement with previous literature data. For the basal plane, as well as for the prism plane, no evident dependence on the ordering of water molecules that constitute the ice block on the rate of crystallization can be observed.

Enhanced Sampling of Phase Transitions in Coarse-Grained Lipid Bilayers

Freezing and melting of dipalmitoylphosphatidylcholine (DPPC) bilayers are simulated in both the explicit (Wet) and implicit solvent (Dry) coarse-grained MARTINI force fields with enhanced sampling, via the isobaric, molecular dynamics version of the generalized replica exchange method (gREM). Phase transitions are described with the entropic viewpoint, based upon the statistical temperature as a function of enthalpy, T_S(H) = 1/(dS(H)/dH), where S is the configurational entropy. Bilayer thickness, area per lipid, and the second-rank order parameter (P2) are calculated vs temperature in the transition range. In a 32-lipid Wet MARTINI system, transitions in the lipid and water subsystems are strongly coupled, giving rise to considerable structure in T_S(H) and the need to specify the state of the water when reporting a lipid transition temperature. For gel lipid + liquid water → fluid lipid + liquid water, we find 292.4 K. The small system is influenced by finite-size effects, but it is argued that the entropic approach is well suited to revealing them, which will be particularly relevant for studies of finite nanosystems where there is no thermodynamic limit. In a 390-lipid Dry MARTINI system, two-dimensional analogues of the topographies of coexisting states ("subphases") seen in pure fluids are found. They are not seen in the 32-lipid Wet or Dry system, but the Dry lipids show a new type of state with gel in one leaflet and tilted gel in the other. Dry bilayer transition temperatures are 333.3 K (390 lipids) and 338 K (32 lipids), indicating that the 32-lipid system is not too small for a qualitative study of the transition. Physical arguments are given for Dry lipid system size dependence and for the difference between Wet and Dry systems.

Ice Nucleation Efficiency of Hydroxylated Organic Surfaces Is Controlled by Their Structural Fluctuations and Mismatch to Ice

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.

Free energy contributions and structural characterization of stacking disordered ices

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.