Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates

Dissociation line and driving force for nucleation of the nitrogen hydrate from computer simulation. II. Effect of multiple occupancy

In this work, we determine the dissociation line of the nitrogen (N2) hydrate by computer simulation using the TIP4P/Ice model for water and the TraPPE force field for N2. This work is the natural extension of Paper I, in which the dissociation temperature of the N2 hydrate has been obtained at 500, 1000, and 1500 bar [Algaba et al., J. Chem. Phys. 159, 224707 (2023)] using the solubility method and assuming single occupancy. We extend our previous study and determine the dissociation temperature of the N2 hydrate at different pressures, from 500 to 4500 bar, taking into account the single and double occupancy of the N2 molecules in the hydrate structure. We calculate the solubility of N2 in the aqueous solution as a function of temperature when it is in contact with a N2-rich liquid phase and when in contact with the hydrate phase with single and double occupancy via planar interfaces. Both curves intersect at a certain temperature that determines the dissociation temperature at a given pressure. We observe a negligible effect of occupancy on the dissociation temperature. Our findings are in very good agreement with the experimental data taken from the literature. We have also obtained the driving force for the nucleation of the hydrate as a function of temperature and occupancy at several pressures. As in the case of the dissociation line, the effect of occupancy on the driving force for nucleation is negligible. To the best of our knowledge, this is the first time that the effect of the occupancy on the driving force for nucleation of a hydrate that exhibits sII crystallographic structure is studied from computer simulation.

Competing Reaction Mechanisms of Peptide Bond Formation in Water Revealed by Deep Potential Molecular Dynamics and Path Sampling

The formation of an amide bond is an essential step in the synthesis of materials and drugs, and in the assembly of amino acids to form peptides. The mechanism of this reaction has been studied extensively, in particular to understand how it can be catalyzed, but a representation capable of explaining all the experimental data is still lacking. Numerical simulation should provide the necessary molecular description, but the solvent involvement poses a number of challenges. Here, we combine the efficiency and accuracy of neural network potential-based reactive molecular dynamics with the extensive and unbiased exploration of reaction pathways provided by transition path sampling. Using microsecond-scale simulations at the density functional theory level, we show that this method reveals the presence of two competing distinct mechanisms for peptide bond formation between alanine esters in aqueous solution. We describe how both reaction pathways, via a general base catalysis mechanism and via direct cleavage of the tetrahedral intermediate respectively, change with pH. This result contrasts with the conventional mechanism involving a single pathway in which only the barrier heights are affected by pH. We show that this new proposal involving two competing mechanisms is consistent with the experimental data, and we discuss the implications for peptide bond formation under prebiotic conditions and in the ribosome. Our work shows that integrating deep potential molecular dynamics with path sampling provides a powerful approach for exploring complex chemical mechanisms.

Enhanced Sampling of Configuration and Path Space in a Generalized Ensemble by Shooting Point Exchange

The computer simulation of many molecular processes is complicated by long timescales caused by rare transitions between long-lived states. Here, we propose a new approach to simulate such rare events, which combines transition path sampling with enhanced exploration of configuration space. The method relies on exchange moves between configuration and trajectory space, carried out based on a generalized ensemble. This scheme substantially enhances the efficiency of the transition path sampling simulations, particularly for systems with multiple transition channels, and yields information on thermodynamics, kinetics and reaction coordinates of molecular processes without distorting their dynamics. The method is illustrated using the isomerization of proline in the KPTP tetrapeptide.

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.

Probing the critical nucleus size in tetrahydrofuran clathrate hydrate formation using surface-anchored nanoparticles

Controlling the formation of clathrate hydrates is crucial for advancing hydrate-based technologies. However, the microscopic mechanism underlying clathrate hydrate formation through nucleation remains poorly elucidated. Specifically, the critical nucleus, theorized as a pivotal step in nucleation, lacks empirical validation. Here, we report uniform nanoparticles, e.g., graphene oxide (GO) nanosheets and gold or silver nanocubes with controlled sizes, as nanoprobes to experimentally measure the size of the critical nucleus of tetrahydrofuran (THF) clathrate hydrate formation. The capability of the nanoparticles in facilitating THF clathrate hydrate nucleation displays generally an abrupt change at a nanoparticle-size-determined specific supercooling. It is revealed that the free-energy barrier shows an abrupt change when the nanoparticles have an approximately the same size as that of the critical nucleus. Thus, it is inferred that THF clathrate hydrate nucleation involves the creation of a critical nucleus with its size being inversely proportional to the supercooling. By proving the existence and determining the supercooling-dependent size of the critical nucleus of the THF clathrate hydrates, this work brings insights in the microscopic pictures of the clathrate hydrate nucleation.

Heterogeneous nucleation in the random field Ising model

We investigate the nucleation dynamics of the three-dimensional random field Ising model under an external field. We use umbrella sampling to compute the free-energy cost of a critical nucleus and use forward flux sampling for the direct estimation of nucleation rates. For moderate to strong disorder, our results indicate that the size of the nucleating cluster is not a good reaction coordinate, contrary to the pure Ising model. We rectify this problem by introducing a coordinate that also accounts for the location of the nucleus. Using the free energy barrier to predict the nucleation rate, we find reasonable agreement, although deviations become stronger as disorder increases. We attribute this effect to cluster shape fluctuations. We also discuss finite-size effects on the nucleation rate.

Resolving Temperature-Dependent Hydrate Nucleation Pathway: The Role of "Transition Layer"

Understanding the nucleation of natural gas hydrate (NGH) at different conditions has important implications to NGH recovery and other industrial applications, such as gas storage and separation. Herein, vast numbers of hydrate nucleation events are traced via molecular dynamics (MD) simulations at different degrees of supercooling (or driving forces). Specifically, to precisely characterize a hydrate nucleus from an aqueous system during the MD simulation, we develop an evolutionary order parameter (OP) to recognize the nucleus size and shape. Subsequently, the free energy landscapes of hydrate during nucleation are explored by using the newly developed OP. The results suggest that at 270 K (or 0.92 Tm supercooling, where Tm is the melting point), the near-rounded nucleus prevails during the nucleation, as described from the classical nucleation theory. In contrast, at relatively strong driving forces of 0.85 and 0.88 Tm, nonclassical nucleation events arise. Specifically, the pathway toward an elongated nucleus becomes as important as the pathway toward a near-rounded nucleus. To explain the distinct nucleation phenomena at different supercoolings, a notion of a "transition layer" (or liquid-blob-like layer) is proposed. Here, the transition layer is to describe the interfacial region between the nucleus and aqueous solution, and this layer entails two functionalities: (1) it tends to retain CH4 depending on the degrees of supercooling and (2) it facilitates collision among CH4, which thus promote the incorporation of CH4 into nucleus. Our simulation indicates that compared to the near-rounded nucleus, the transition layer surrounding the elongated nucleus is more evident with the higher collision rate among CH4 molecules. As such, the transition layer tends to promote the elongated nucleus pathway, while offsetting the cost of larger surface free energy associated with the elongated nucleus. At 0.92 Tm, however, the transition layer gradually disappears, and classical nucleation events dominate. Overall, the notion of "transition layer" offers deeper insight into the NGH nucleation at different degrees of supercooling and could be extended to describe other types of hydrate nucleation.

Microscopic insights on clathrate hydrate growth from non-equilibrium molecular dynamics simulations

Clathrate hydrates form and grow at interfaces. Understanding the relevant molecular processes is crucial for developing hydrate-based technologies. Many computational studies focus on hydrate growth within the aqueous phase using the 'direct coexistence method', which is limited in its ability to investigate hydrate film growth at hydrocarbon-water interfaces. To overcome this shortcoming, a new simulation setup is presented here, which allows us to study the growth of a methane hydrate nucleus in a system where oil-water, hydrate-water, and hydrate-oil interfaces are all simultaneously present, thereby mimicking experimental setups. Using this setup, hydrate growth is studied here under the influence of two additives, a polyvinylcaprolactam oligomer and sodium dodecyl sulfate, at varying concentrations. Our results confirm that hydrate films grow along the oil-water interface, in general agreement with visual experimental observations; growth, albeit slower, also occurs at the hydrate-water interface, the interface most often interrogated via simulations. The results obtained demonstrate that the additives present within curved interfaces control the solubility of methane in the aqueous phase, which correlates with hydrate growth rate. Building on our simulation insights, we suggest that by combining data for the potential of mean force profile for methane transport across the oil-water interface and for the average free energy required to perturb a flat interface, it is possible to predict the performance of additives used to control hydrate growth. These insights could be helpful to achieve optimal methane storage in hydrates, one of many applications which are attracting significant fundamental and applied interests.

Solubility of carbon dioxide in water: Some useful results for hydrate nucleation

In this paper, the solubility of carbon dioxide (CO2) in water along the isobar of 400 bar is determined by computer simulations using the well-known TIP4P/Ice force field for water and the TraPPE model for CO2. In particular, the solubility of CO2 in water when in contact with the CO2 liquid phase and the solubility of CO2 in water when in contact with the hydrate have been determined. The solubility of CO2 in a liquid-liquid system decreases as the temperature increases. The solubility of CO2 in a hydrate-liquid system increases with temperature. The two curves intersect at a certain temperature that determines the dissociation temperature of the hydrate at 400 bar (T3). We compare the predictions with T3 obtained using the direct coexistence technique in a previous work. The results of both methods agree, and we suggest 290(2) K as the value of T3 for this system using the same cutoff distance for dispersive interactions. We also propose a novel and alternative route to evaluate the change in chemical potential for the formation of hydrates along the isobar. The new approach is based on the use of the solubility curve of CO2 when the aqueous solution is in contact with the hydrate phase. It considers rigorously the non-ideality of the aqueous solution of CO2, providing reliable values for the driving force for nucleation of hydrates in good agreement with other thermodynamic routes used. It is shown that the driving force for hydrate nucleation at 400 bar is larger for the methane hydrate than for the carbon dioxide hydrate when compared at the same supercooling. We have also analyzed and discussed the effect of the cutoff distance of dispersive interactions and the occupancy of CO2 on the driving force for nucleation of the hydrate.

Machine-guided path sampling to discover mechanisms of molecular self-organization

Molecular self-organization driven by concerted many-body interactions produces the ordered structures that define both inanimate and living matter. Here we present an autonomous path sampling algorithm that integrates deep learning and transition path theory to discover the mechanism of molecular self-organization phenomena. The algorithm uses the outcome of newly initiated trajectories to construct, validate and-if needed-update quantitative mechanistic models. Closing the learning cycle, the models guide the sampling to enhance the sampling of rare assembly events. Symbolic regression condenses the learned mechanism into a human-interpretable form in terms of relevant physical observables. Applied to ion association in solution, gas-hydrate crystal formation, polymer folding and membrane-protein assembly, we capture the many-body solvent motions governing the assembly process, identify the variables of classical nucleation theory, uncover the folding mechanism at different levels of resolution and reveal competing assembly pathways. The mechanistic descriptions are transferable across thermodynamic states and chemical space.

Molecular Dynamics Study of Clathrate-like Ordering of Water in Supersaturated Methane Solution at Low Pressure

Using molecular dynamics, the evolution of a metastable solution for “methane + water” was studied for concentrations of 3.36, 6.5, 9.45, 12.2, and 14.8 mol% methane at 270 K and 1 bar during 100 ns. We have found the intriguing behavior of the system containing over 10,000 water molecules: the formation of hydrate-like structures is observed at 6.5 and 9.45 mol% concentrations throughout the entire solution volume. This formation of “blobs” and the following amorphous hydrate were studied. The creation of a metastable methane solution through supersaturation is the key to triggering the collective process of hydrate formation under low pressure. Even the first stage (0–1 ns), before the first fluctuating cavities appear, is a collective process of H-bond network reorganization. The formation of fluctuation cavities appears before steady hydrate growth begins and is associated with a preceding uniform increase in the water molecule’s tetrahedrality. Later, the constantly presented hydrate cavities become the foundation for a few independent hydrate nucleation centers, this evolution is consistent with the labile cluster and local structure hypotheses. This new mechanism of hydrogen-bond network reorganization depends on the entropy of the cavity arrangement of the guest molecules in the hydrate lattice and leads to hydrate growth.

Homogeneous nucleation rate of methane hydrate formation under experimental conditions from seeding simulations

In this work, we shall estimate via computer simulations the homogeneous nucleation rate for the methane hydrate at 400 bars for a supercooling of about 35 K. The TIP4P/ICE model and a Lennard-Jones center were used for water and methane, respectively. To estimate the nucleation rate, the seeding technique was employed. Clusters of the methane hydrate of different sizes were inserted into the aqueous phase of a two-phase gas-liquid equilibrium system at 260 K and 400 bars. Using these systems, we determined the size at which the cluster of the hydrate is critical (i.e., it has 50% probability of either growing or melting). Since nucleation rates estimated from the seeding technique are sensitive to the choice of the order parameter used to determine the size of the cluster of the solid, we considered several possibilities. We performed brute force simulations of an aqueous solution of methane in water in which the concentration of methane was several times higher than the equilibrium concentration (i.e., the solution was supersaturated). From brute force runs, we infer the value of the nucleation rate for this system rigorously. Subsequently, seeding runs were carried out for this system, and it was found that only two of the considered order parameters were able to reproduce the value of the nucleation rate obtained from brute force simulations. By using these two order parameters, we estimated the nucleation rate under experimental conditions (400 bars and 260 K) to be of the order of log₁₀ (J/(m³ s)) = -7(5).

Homogeneous nucleation of crystalline methane hydrate in molecular dynamics transition paths sampled under realistic conditions

Methane hydrates are important from a scientific and industrial perspective, and form by nucleation and growth from a supersaturated aqueous solution of methane. Molecular simulation is able to shed light on the process of homogeneous nucleation of hydrates, using straightforward molecular dynamics or rare event enhanced sampling techniques with atomistic and coarse grained force fields. In our previous work [Arjun, T. A. Berendsen, and P. G. Bolhuis, Proc. Natl. Acad. Sci. U. S. A. 116, 19305 (2019)], we performed transition path sampling (TPS) simulations using all atom force fields under moderate driving forces at high pressure, which enabled unbiased atomistic insight into the formation of methane hydrates. The supersaturation in these simulations was influenced by the Laplace pressure induced by the spherical gas reservoir. Here, we investigate the effect of removing this influence. Focusing on the supercooled, supersaturated regime to keep the system size tractable, our TPS simulations indicate that nuclei form amorphous structures below roughly 260 K and crystalline sI structures above 260 K. For these temperatures, the average transition path lengths are significantly longer than in our previous study, pushing the boundaries of what can be achieved with TPS. The temperature to observe a critical nucleus of certain size was roughly 20 K lower compared to a spherical reservoir due to the lower concentration of methane in the solution, yielding a reduced driving force. We analyze the TPS results using a model based on classical nucleation theory. The corresponding free energy barriers are estimated and found to be consistent with previous predictions, thus adding to the overall picture of the hydrate formation process.

Critical comparison of general-purpose collective variables for crystal nucleation

The nucleation of crystals is a prominent phenomenon in science and technology that still lacks a full atomic-scale understanding. Much work has been devoted to identifying order parameters able to track the process, from the inception of early nuclei to their maturing to critical size until growth of an extended crystal. We critically assess and compare two powerful distance-based collective variables, an effective entropy derived from liquid state theory and the path variable based on permutation invariant vectors using the Kob-Andersen binary mixture and a combination of enhanced-sampling techniques. Our findings reveal a comparable ability to drive nucleation when a bias potential is applied, and comparable free-energy barriers and structural features. Yet, we also found an imperfect correlation with the committor probability on the barrier top which was bypassed by changing the order parameter definition.

Effects of marine environments on methane hydrate formation in clay nanopores: A molecular dynamics study

In nature, CH₄ hydrates are mainly buried in marine sediments. The complex marine environments on the seafloor continuously affect hydrate formation. Herein, systematic molecular simulations have been performed to investigate CH₄ hydrate formation in clay nanopore, mainly affected by several marine environmental factors, including seawater salinity, pressure and temperature. Simulation results show that these factors exert different effects on hydrate formation in the nanopore and the outside bulk solutions by affecting the mass transfer and phase separation inside and outside of the nanopore. Specifically, high salinity hinders the diffusion of CH₄ molecules from nanopores into the outside bulk solutions, promoting hydrate formation in nanopore and inhibiting hydrate formation in bulk solution; salinity has a dual effect on hydrate formation in the whole system by changing the local CH₄ concentration via the formation of the hydration of salt ions. High pressure favors the diffusion of CH₄ molecules from nanopore into outside bulk solutions, promoting hydrate formation in bulk solution and inhibiting hydrate formation in nanopore; high pressure promotes hydrate formation at the nanopore throats by increasing CH₄ concentration and reducing ion concentration therein. In contrast, temperature significantly affects hydrate formation in the system by causing phase separation, i.e. high temperature promotes the aggregation of CH₄ molecules to form nanobubbles and inhibits hydrate formation. Under high temperature conditions, the nanobubble in the nanopore gradually decomposes, while the nanobubble in the outside bulk solution grows an extra-large cylindrical nanobubble. These molecular insights into the formation behavior of CH₄ hydrates in clay nanopores are helpful for understanding the formation process of natural gas hydrates in marine sediments and the development and utilization of CH₄ hydrates.

Solubility of Methane in Water: Some Useful Results for Hydrate Nucleation

In this paper, the solubility of methane in water along the 400 bar isobar is determined by computer simulations using the TIP4P/Ice force field for water and a simple LJ model for methane. In particular, the solubility of methane in water when in contact with the gas phase and the solubility of methane in water when in contact with the hydrate has been determined. The solubility of methane in a gas-liquid system decreases as temperature increases. The solubility of methane in a hydrate-liquid system increases with temperature. The two curves intersect at a certain temperature that determines the triple point T₃ at a certain pressure. We also determined T₃ by the three-phase direct coexistence method. The results of both methods agree, and we suggest 295(2) K as the value of T₃ for this system. We also analyzed the impact of curvature on the solubility of methane in water. We found that the presence of curvature increases the solubility in both the gas-liquid and hydrate-liquid systems. The change in chemical potential for the formation of hydrate is evaluated along the isobar using two different thermodynamic routes, obtaining good agreement between them. It is shown that the driving force for hydrate nucleation under experimental conditions is higher than that for the formation of pure ice when compared at the same supercooling. We also show that supersaturation (i.e., concentrations above those of the planar interface) increases the driving force for nucleation dramatically. The effect of bubbles can be equivalent to that of an additional supercooling of about 20 K. Having highly supersaturated homogeneous solutions makes possible the spontaneous formation of the hydrate at temperatures as high as 285 K (i.e., 10K below T₃). The crucial role of the concentration of methane for hydrate formation is clearly revealed. Nucleation of the hydrate can be either impossible or easy and fast depending on the concentration of methane which seems to play the leading role in the understanding of the kinetics of hydrate formation.

Mixed Hydrate Nucleation: Molecular Mechanisms and Cage Structures

The molecular-level details of the formation of mixed gas hydrates remain elusive despite their significance for a variety of scientific and industrial applications. In this study, extensive molecular simulations have been performed to examine the behavior of CH₄/H₂S mixed hydrate nucleation utilizing two different simulation setups varying in compositions and temperatures. The observed behavior exhibits similar phenomenology across the various systems once differences in nucleation rates and guest uptake are accounted for. We find that CH₄ is always enriched in the hydrate phase while the aqueous phase is enriched in H₂S. Even with H₂S as a minor component (i.e., 10% mole fraction), the system can mirror the overall nucleation kinetics of pure H₂S hydrate systems with CH₄-dominant nuclei. Through analyses of cages and their transitions, nonstandard cages, particularly those with 12 faces (e.g., 5¹⁰6²), have been found to be key intermediate cage types in the early stage of nucleation. Additionally, we present previously unreported cage types comprising heptagonal faces (e.g., 5⁹6²7¹) as having a significant role in the early-stage gas hydrate structural transitions.

Double Life of Methanol: Experimental Studies and Nonequilibrium Molecular-Dynamics Simulation of Methanol Effects on Methane-Hydrate Nucleation

We have investigated systematically and statistically methanol-concentration effects on methane-hydrate nucleation using both experiment and restrained molecular-dynamics simulation, employing simple observables to achieve an initially homogeneous methane-supersaturated solution particularly favorable for nucleation realization in reasonable simulation times. We observe the pronounced "bifurcated" character of the nucleation rate upon methanol concentration in both experiments and simulation, with promotion at low concentrations and switching to industrially familiar inhibition at higher concentrations. Higher methanol concentrations suppress hydrate growth by in-lattice methanol incorporation, resulting in the formation of "defects", increasing the energy of the nucleus. At low concentrations, on the contrary, the detrimental effect of defects is more than compensated for by the beneficial contribution of CH₃ in easing methane incorporation in the cages or replacing it altogether.

Explaining reaction coordinates of alanine dipeptide isomerization obtained from deep neural networks using Explainable Artificial Intelligence (XAI)

A method for obtaining appropriate reaction coordinates is required to identify transition states distinguishing product and reactant in complex molecular systems. Recently, abundant research has been devoted to obtaining reaction coordinates using artificial neural networks from deep learning literature, where many collective variables are typically utilized in the input layer. However, it is difficult to explain the details of which collective variables contribute to the predicted reaction coordinates owing to the complexity of the nonlinear functions in deep neural networks. To overcome this limitation, we used Explainable Artificial Intelligence (XAI) methods of the Local Interpretable Model-agnostic Explanation (LIME) and the game theory-based framework known as Shapley Additive exPlanations (SHAP). We demonstrated that XAI enables us to obtain the degree of contribution of each collective variable to reaction coordinates that is determined by nonlinear regressions with deep learning for the committor of the alanine dipeptide isomerization in vacuum. In particular, both LIME and SHAP provide important features to the predicted reaction coordinates, which are characterized by appropriate dihedral angles consistent with those previously reported from the committor test analysis. The present study offers an AI-aided framework to explain the appropriate reaction coordinates, which acquires considerable significance when the number of degrees of freedom increases.

An extended autoencoder model for reaction coordinate discovery in rare event molecular dynamics datasets

The reaction coordinate (RC) is the principal collective variable or feature that determines the progress along an activated or reactive process. In a molecular simulation using enhanced sampling, a good description of the RC is crucial for generating sufficient statistics. Moreover, the RC provides invaluable atomistic insight into the process under study. The optimal RC is the committor, which represents the likelihood of a system to evolve toward a given state based on the coordinates of all its particles. As the interpretability of such a high dimensional function is low, a more practical approach is to describe the RC by some low-dimensional molecular collective variables or order parameters. While several methods can perform this dimensionality reduction, they usually require a preselection of these low-dimension collective variables (CVs). Here, we propose to automate this dimensionality reduction using an extended autoencoder, which maps the input (many CVs) onto a lower-dimensional latent space, which is subsequently used for the reconstruction of the input as well as the prediction of the committor function. As a consequence, the latent space is optimized for both reconstruction and committor prediction and is likely to yield the best non-linear low-dimensional representation of the committor. We test our extended autoencoder model on simple but nontrivial toy systems, as well as extensive molecular simulation data of methane hydrate nucleation. The extended autoencoder model can effectively extract the underlying mechanism of a reaction, make reliable predictions about the committor of a given configuration, and potentially even generate new paths representative for a reaction.

The seven deadly sins: When computing crystal nucleation rates, the devil is in the details

The formation of crystals has proven to be one of the most challenging phase transformations to quantitatively model-let alone to actually understand-be it by means of the latest experimental technique or the full arsenal of enhanced sampling approaches at our disposal. One of the most crucial quantities involved with the crystallization process is the nucleation rate, a single elusive number that is supposed to quantify the average probability for a nucleus of critical size to occur within a certain volume and time span. A substantial amount of effort has been devoted to attempt a connection between the crystal nucleation rates computed by means of atomistic simulations and their experimentally measured counterparts. Sadly, this endeavor almost invariably fails to some extent, with the venerable classical nucleation theory typically blamed as the main culprit. Here, we review some of the recent advances in the field, focusing on a number of perhaps more subtle details that are sometimes overlooked when computing nucleation rates. We believe it is important for the community to be aware of the full impact of aspects, such as finite size effects and slow dynamics, that often introduce inconspicuous and yet non-negligible sources of uncertainty into our simulations. In fact, it is key to obtain robust and reproducible trends to be leveraged so as to shed new light on the kinetics of a process, that of crystal nucleation, which is involved into countless practical applications, from the formulation of pharmaceutical drugs to the manufacturing of nano-electronic devices.

Free energy surface of two-step nucleation

We test the theoretical free energy surface (FES) for two-step nucleation (TSN) proposed by Iwamatsu [J. Chem. Phys. 134, 164508 (2011)] by comparing the predictions of the theory to numerical results for the FES recently reported from Monte Carlo simulations of TSN in a simple lattice system [James et al., J. Chem. Phys. 150, 074501 (2019)]. No adjustable parameters are used to make this comparison. That is, all the parameters of the theory are evaluated directly for the model system, yielding a predicted FES, which we then compare to the FES obtained from simulations. We find that the theoretical FES successfully predicts the numerically evaluated FES over a range of thermodynamic conditions that spans distinct regimes of behavior associated with TSN. All the qualitative features of the FES are captured by the theory, and the quantitative comparison is also very good. Our results demonstrate that Iwamatsu's extension of classical nucleation theory provides an excellent framework for understanding the thermodynamics of TSN.

Homogenous nucleation rate of CO2 hydrates using transition interface sampling

Carbon dioxide and water can form solid clathrate structures in which water cages encapsulate the gas molecules. Such hydrates have sparked much interest due to their possible application in CO₂ sequestration. How the solid structure forms exactly from the liquid phase via a homogenous nucleation process is still poorly understood. This nucleation event is rare on the molecular timescale even under moderate undercooling or supersaturation conditions because of the large free energy barrier toward crystallization, rendering a brute force simulation of hydrate nucleation unfeasible for moderate undercooling or supersaturation. Here, we perform transition interface sampling simulations to quantify the homogenous nucleation rate for CO₂ hydrate formation using accurate atomistic force fields at 500 bars for three different temperatures between 260 and 273 K. Collecting more than 100 000 pathways comprising roughly two milliseconds of simulation time, we computed a nucleation rate in the amorphous phase of ∼10²¹ nuclei s^-1 cm^-3 for a temperature of 260 K and a rate of ∼10¹² nuclei s^-1 cm^-3 for a temperature of 265 K. For a temperature of 273 K, we find that the hydrate forms an sI crystalline phase with a rate of order of ∼10¹ nuclei s^-1 cm^-3. We compare these rates to classical nucleation theory estimates as well as experiments, and to nucleation rate estimates for methane hydrates and discuss possible causes of the observed differences. Our findings shed light on the kinetics of this important clathrate and should assist in future hydrate formation investigation.

Reduced phase stability and faster formation/dissociation kinetics in confined methane hydrate

The mechanisms involved in the formation/dissociation of methane hydrate confined at the nanometer scale are unraveled using advanced molecular modeling techniques combined with a mesoscale thermodynamic approach. Using atom-scale simulations probing coexistence upon confinement and free energy calculations, phase stability of confined methane hydrate is shown to be restricted to a narrower temperature and pressure domain than its bulk counterpart. The melting point depression at a given pressure, which is consistent with available experimental data, is shown to be quantitatively described using the Gibbs-Thomson formalism if used with accurate estimates for the pore/liquid and pore/hydrate interfacial tensions. The metastability barrier upon hydrate formation and dissociation is found to decrease upon confinement, therefore providing a molecular-scale picture for the faster kinetics observed in experiments on confined gas hydrates. By considering different formation mechanisms-bulk homogeneous nucleation, external surface nucleation, and confined nucleation within the porosity-we identify a cross-over in the nucleation process; the critical nucleus formed in the pore corresponds either to a hemispherical cap or to a bridge nucleus depending on temperature, contact angle, and pore size. Using the classical nucleation theory, for both mechanisms, the typical induction time is shown to scale with the pore volume to surface ratio and hence the pore size. These findings for the critical nucleus and nucleation rate associated with such complex transitions provide a means to rationalize and predict methane hydrate formation in any porous media from simple thermodynamic data.

Parasitic crystallization of colloidal electrolytes: growing a metastable crystal from the nucleus of a stable phase

Colloidal particles have been extensively used to comprehend the main principles governing liquid-crystal nucleation. Multiple mechanisms and frameworks have been proposed, through either experiments or computational approaches, to rationalise the ubiquitous formation of colloidal crystals. In this work, we elucidate the nucleation scenario behind the crystallization of oppositely charged colloids. By performing molecular dynamics simulations of colloidal electrolytes in combination with the Seeding technique, we evaluate the fundamental factors, such as the nucleation rate, free energy barrier, surface tension and kinetic pre-factor, that determine the liquid-to-solid transition of several crystalline polymorphs. Our results show that at a high packing fraction, there is a cross-over between the nucleation of the CsCl structure and that of a substitutionally disordered fcc phase, despite the CuAu crystal being the most stable phase. We demonstrate that the crucial factor in determining which phase nucleates the fastest is the free energy cost of the cluster formation rather than their kinetic ability to grow from the liquid. While at a low packing fraction, the stable phase, CsCl, is the one that nucleates and subsequently grows, we show how at moderate and high packing fractions, a disordered fcc phase subsequently grows regardless of the nature of the nucleating phase, termed parasitic crystallization. Taken together, our results provide a panoramic perspective of the complex nucleation scenario of oppositely charged colloids at moderate temperature and rationalise the different thermodynamic and kinetic aspects behind it.

Molecular Understanding of Homogeneous Nucleation of CO2 Hydrates Using Transition Path Sampling

Carbon dioxide hydrate is a solid built from hydrogen-bond stabilized water cages that encapsulate individual CO2 molecules. As potential candidates for reducing greenhouse gases, hydrates have attracted attention from both the industry and scientific community. Under high pressure and low temperature, hydrates are formed spontaneously from a mixture of CO2 and water via nucleation and growth. Yet, for moderate undercooling, i.e., moderate supersaturation, studying hydrate formation with molecular simulations is very challenging due to the high nucleation barriers involved. We investigate the homogeneous nucleation mechanism of CO2 hydrate as a function of temperature using transition path sampling (TPS), which generates ensembles of unbiased dynamical trajectories across the high barrier between the liquid and solid states. The resulting path ensembles reveal that at high driving force (low temperature), amorphous structures are predominantly formed, with 4151062 cages being the most abundant. With increasing temperature, the nucleation mechanism changes, and 51262 becomes the most abundant cage type, giving rise to the crystalline sI structure. Reaction coordinate analysis can reveal the most important collective variable involved in the mechanism. With increasing temperature, we observe a shift from a single feature (size of the nucleus) to a 2-dimensional (size and cage type) variable as the salient ingredient of the reaction coordinate, and then back to only the nucleus size. This finding is in line with the underlying shift from an amorphous to a crystalline nucleation channel. Modeling such complex phase transformations using transition path sampling gives unbiased insight into the molecular mechanisms toward different polymorphs, and how these are determined by thermodynamics and kinetics. This study will be beneficial for researchers aiming to produce such hydrates with different polymorphic forms.

Rate Prediction for Homogeneous Nucleation of Methane Hydrate at Moderate Supersaturation Using Transition Interface Sampling

The crystallization of methane hydrates via homogeneous nucleation under natural, moderate conditions is of both industrial and scientific relevance, yet still poorly understood. Predicting the nucleation rates at such conditions is notoriously difficult due to high nucleation barriers, and requires, besides an accurate molecular model, enhanced sampling. Here, we apply the transition interface sampling technique, which efficiently computes the exact rate of nucleation by generating ensembles of unbiased dynamical trajectories crossing predefined interfaces located between the stable states. Using an accurate atomistic force field and focusing on specific conditions of 280 K and 500 bar, we compute for nucleation directly into the sI crystal phase at a rate of ∼10^-17 nuclei per nanosecond per simulation volume or ∼10² nuclei per second per cm³, in agreement with consensus estimates for nearby conditions. As this is most likely fortuitous, we discuss the causes of the large differences between our results and previous simulation studies. Our work shows that it is now possible to compute rates for methane hydrates at moderate supersaturation, without relying on any assumptions other than the force field.

Propane and Water: The Cooperativity of Unlikely Molecules to Form Clathrate Structures

Many unanswered questions still exist at the molecular level to understand the nucleation process and mechanism of clathrate hydrates, especially for larger guest molecules that would result in the structure II crystal. Here, we report on molecular dynamics simulations for propane and water to describe the molecular mechanism leading to a structure II system. Through a large number (30) of long (5 μs) and coupled annealing (20 μs) simulations, we detail the prenucleation, nucleation, growth, and annealing of propane clathrate hydrate structures at 250 K and 1800 bar. The results demonstrate the equal importance of the empty and occupied cages in the nucleation of propane hydrates. The critical nucleus size is identified to be eight cages. While separate distinct clusters may exist during the prenucleation period, only one survives to grow beyond the critical nucleus size, with the others remaining subcritical. From the annealing simulations, it is clear that solid rearrangement is a very slow process, and 20 μs is still not long enough to capture long-range ordering resembling the structure II crystal. These results, along with the developed analysis method, have a significant impact in advancing our understanding of the nucleation process for unlike molecules.