Optimizing Nucleus Size Metrics for Liquid-Solid Nucleation from Transition Paths of Near-Nanosecond Duration

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.

Divergence among Local Structure, Dynamics, and Nucleation Outcome in Heterogeneous Nucleation of Close-Packed Crystals

Heterogeneous crystal nucleation is the dominant mechanism of crystallization in most systems, yet its underlying physics remains an enigma. While emergent interfacial crystalline order precedes heterogeneous nucleation, its importance in the nucleation mechanism is unclear. Here, we use path sampling simulations of two model systems to demonstrate that crystalline order in its traditional sense is not predictive of the outcome of the heterogeneous nucleation of close-packed crystals. Consequently, structure-based collective variables (CVs) that reliably describe homogeneous nucleation can be poor descriptors of heterogeneous nucleation. This divergence between structure and nucleation outcome is accompanied by an intriguing dynamical anomaly, wherein low-coordinated crystalline particles outpace their liquid-like counterparts. We use committor analysis, high-throughput screening, and machine learning to devise CV optimization strategies and present suitable structural heuristics within the metastable fluid for CV prescreening. Employing such optimized CVs is pivotal for properly characterizing the mechanism of heterogeneous nucleation in metallic and colloidal systems.

An exploration of machine learning models for the determination of reaction coordinates associated with conformational transitions

Determining collective variables (CVs) for conformational transitions is crucial to understanding their dynamics and targeting them in enhanced sampling simulations. Often, CVs are proposed based on intuition or prior knowledge of a system. However, the problem of systematically determining a proper reaction coordinate (RC) for a specific process in terms of a set of putative CVs can be achieved using committor analysis (CA). Identifying essential degrees of freedom that govern such transitions using CA remains elusive because of the high dimensionality of the conformational space. Various schemes exist to leverage the power of machine learning (ML) to extract an RC from CA. Here, we extend these studies and compare the ability of 17 different ML schemes to identify accurate RCs associated with conformational transitions. We tested these methods on an alanine dipeptide in vacuum and on a sarcosine dipeptoid in an implicit solvent. Our comparison revealed that the light gradient boosting machine method outperforms other methods. In order to extract key features from the models, we employed Shapley Additive exPlanations analysis and compared its interpretation with the "feature importance" approach. For the alanine dipeptide, our methodology identifies ϕ and θ dihedrals as essential degrees of freedom in the C7ax to C7eq transition. For the sarcosine dipeptoid system, the dihedrals ψ and ω are the most important for the cisαD to transαD transition. We further argue that analysis of the full dynamical pathway, and not just endpoint states, is essential for identifying key degrees of freedom governing transitions.

A variational approach to assess reaction coordinates for two-step crystallization

Molecule- and particle-based simulations provide the tools to test, in microscopic detail, the validity of classical nucleation theory. In this endeavor, determining nucleation mechanisms and rates for phase separation requires an appropriately defined reaction coordinate to describe the transformation of an out-of-equilibrium parent phase for which myriad options are available to the simulator. In this article, we describe the application of the variational approach to Markov processes to quantify the suitability of reaction coordinates to study crystallization from supersaturated colloid suspensions. Our analysis indicates that collective variables (CVs) that correlate with the number of particles in the condensed phase, the system potential energy, and approximate configurational entropy often feature as the most appropriate order parameters to quantitatively describe the crystallization process. We apply time-lagged independent component analysis to reduce high-dimensional reaction coordinates constructed from these CVs to build Markov State Models (MSMs), which indicate that two barriers separate a supersaturated fluid phase from crystals in the simulated environment. The MSMs provide consistent estimates for crystal nucleation rates, regardless of the dimensionality of the order parameter space adopted; however, the two-step mechanism is only consistently evident from spectral clustering of the MSMs in higher dimensions. As the method is general and easily transferable, the variational approach we adopt could provide a useful framework to study controls for crystal nucleation.

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.

Nonlinear Reaction Coordinate of an Enzyme Catalyzed Proton Transfer Reaction

We present an in-depth study on the theoretical calculation of an optimum reaction coordinate as a linear or nonlinear combination of important collective variables (CVs) sampled from an ensemble of reactive transition paths for an intramolecular proton transfer reaction catalyzed by the enzyme human carbonic anhydrase (HCA) II. The linear models are optimized by likelihood maximization for a given number of CVs. The nonlinear models are based on an artificial neural network with the same number of CVs and optimized by minimizing the root-mean-square error in comparison to a training set of committor estimators generated for the given transition. The nonlinear reaction coordinate thus obtained yields the free energy of activation and rate constant as 9.46 kcal mol^-1 and 1.25 × 10⁶ s^-1, respectively. These estimates are found to be in quantitative agreement with the known experimental results. We have also used an extended autoencoder model to show that a similar analysis can be carried out using a single CV only. The resultant free energies and kinetics of the reaction slightly overestimate the experimental data. The implications of these results are discussed using a detailed microkinetic scheme of the proton transfer reaction catalyzed by HCA II.

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.

Multivalent Surface Cations Enhance Heterogeneous Freezing of Water on Muscovite Mica

Heterogeneous ice nucleation is a crucial phenomenon in various fields of fundamental and applied science. We investigate the effect of surface cations on freezing of water on muscovite mica. Mica is unique in that the exposed ion on its surface can be readily and easily exchanged without affecting other properties such as surface roughness. We investigate freezing on natural (K⁺) mica and mica in which we have exchanged K⁺ for Al³⁺, Mg²⁺, Ca²⁺, and Sr²⁺. We find that liquid water freezes at higher temperatures when ions of higher valency are present on the surface, thus exposing more of the underlying silica layer. Our data also show that the size of the ion affects the characteristic freezing temperature. Using molecular dynamics simulations, we investigate the effects that the ion valency and exposed silica layer have on the behavior of water on the surface. The results indicate that multivalent cations enhance the probability of forming large clusters of hydrogen bonded water molecules that are anchored by the hydration shells of the cations. These clusters also have a large fraction of free water that can reorient to take ice-like configurations, which are promoted by the regions on mica devoid of the ions. Thus, these clusters could serve as seedbeds for ice nuclei. The combined experimental and simulation studies shed new light on the influence of surface ions on heterogeneous ice nucleation.

Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates

Mixtures of methane gas and water can form ice-like solid methane hydrates via homogeneous nucleation at natural, moderate conditions. Understanding the formation of the different hydrate polymorphs, required for improved control of crystallization, is hampered by limited experimental resolution. Direct molecular dynamics simulations could provide insight, but it would require computation times beyond the age of the Universe for a single nucleation event to take place. Yet, the nucleation event itself, while rare, occurs fast. Therefore, we sample and examine ensembles of rare unbiased nucleation trajectories. Detailed analysis shows how selection between competing amorphous and crystalline polymorph formation mechanisms takes place. The conclusions shed light onto the formation of clathrate hydrates.

Methane hydrates have important industrial and climate implications, yet their formation via homogeneous nucleation under natural, moderate conditions is poorly understood. Obtaining such understanding could lead to improved control of crystallization, as well as insight into polymorph selection in general, but is hampered by limited experimental resolution. Direct molecular dynamics simulations using atomistic force fields could provide such insight, but are not feasible for moderate undercooling, due to the rare event nature of nucleation. Instead, we harvest ensembles of the rare unbiased nucleation trajectories by employing transition path sampling. We find that with decreasing undercooling the mechanism shifts from amorphous to crystalline polymorph formation. At intermediate temperature the 2 mechanisms compete. Reaction coordinate analysis reveals the amount of a specific methane cage type is crucial for crystallization, while irrelevant for amorphous solids. Polymorph selection is thus governed by kinetic accessibility of the correct cage type and, moreover, occurs at precritical nucleus sizes, apparently against Ostwald’s step rule. We argue that these results are still in line with classical nucleation theory. Our findings illuminate how selection between competing methane hydrate polymorphs occurs and might generalize to other hydrates and molecular crystal formation.

Molecular Modeling Investigations of Sorption and Diffusion of Small Molecules in Glassy Polymers

With a wide range of applications, from energy and environmental engineering, such as in gas separations and water purification, to biomedical engineering and packaging, glassy polymeric materials remain in the core of novel membrane and state-of the art barrier technologies. This review focuses on molecular simulation methodologies implemented for the study of sorption and diffusion of small molecules in dense glassy polymeric systems. Basic concepts are introduced and systematic methods for the generation of realistic polymer configurations are briefly presented. Challenges related to the long length and time scale phenomena that govern the permeation process in the glassy polymer matrix are described and molecular simulation approaches developed to address the multiscale problem at hand are discussed.

Crystallization of hard spheres revisited. I. Extracting kinetics and free energy landscape from forward flux sampling

We investigate the kinetics and the free energy landscape of the crystallization of hard spheres from a supersaturated metastable liquid though direct simulations and forward flux sampling. In this first paper, we describe and test two different ways to reconstruct the free energy barriers from the sampled steady state probability distribution of cluster sizes without sampling the equilibrium distribution. The first method is based on mean first passage times, and the second method is based on splitting probabilities. We verify both methods for a single particle moving in a double-well potential. For the nucleation of hard spheres, these methods allow us to probe a wide range of supersaturations and to reconstruct the kinetics and the free energy landscape from the same simulation. Results are consistent with the scaling predicted by classical nucleation theory although a quantitative fit requires a rather large effective interfacial tension.

Pathways to self-organization: Crystallization via nucleation and growth

Crystallization, a prototypical self-organization process during which a disordered state spontaneously transforms into a crystal characterized by a regular arrangement of its building blocks, usually proceeds by nucleation and growth. In the initial stages of the transformation, a localized nucleus of the new phase forms in the old one due to a random fluctuation. Most of these nuclei disappear after a short time, but rarely a crystalline embryo may reach a critical size after which further growth becomes thermodynamically favorable and the entire system is converted into the new phase. In this article, we will discuss several theoretical concepts and computational methods to study crystallization. More specifically, we will address the rare event problem arising in the simulation of nucleation processes and explain how to calculate nucleation rates accurately. Particular attention is directed towards discussing statistical tools to analyze crystallization trajectories and identify the transition mechanism.

S-shooting: a Bennett–Chandler-like method for the computation of rate constants from committor trajectories

Mechanisms of rare transitions between long-lived stable states are often analyzed in terms of commitment probabilities, determined from swarms of short molecular dynamics trajectories. Here, we present a computer simulation method to determine rate constants from such short trajectories combined with free energy calculations. The method, akin to the Bennett–Chandler approach for the calculation of reaction rate constants, requires the definition of a valid reaction coordinate and can be applied to both under- and overdamped dynamics. We verify the correctness of the algorithm using a one-dimensional random walker in a double-well potential and demonstrate its applicability to complex transitions in condensed systems by calculating cavitation rates for water at negative pressures.

Molecular simulation of homogeneous nucleation of crystals of an ionic liquid from the melt

The homogeneous nucleation of crystals of the ionic liquid [dmim(+)][Cl(-)] from its supercooled liquid phase in the bulk (P = 1 bar, T = 340 K, representing a supercooling of 58 K) was studied using molecular simulations. The string method in collective variables [Maragliano et al., J. Chem. Phys. 125, 024106 (2006)] was used in combination with Markovian milestoning with Voronoi tessellations [Maragliano et al., J. Chem. Theory Comput. 5, 2589-2594 (2009)] and order parameters for molecular crystals [E. E. Santiso and B. L. Trout, J. Chem. Phys. 134, 064109 (2011)] to sketch a minimum free energy path connecting the supercooled liquid and the monoclinic crystal phases, and to determine the free energy and the rates involved in the homogeneous nucleation process. The physical significance of the configurations found along this minimum free energy path is discussed with the help of calculations based on classical nucleation theory and with additional simulation results obtained for a larger system. Our results indicate that, at a supercooling of 58 K, the liquid has to overcome a free energy barrier of the order of 60 kcal/mol and to form a critical nucleus with an average size of about 3.6 nm, before it reaches the thermodynamically stable crystal phase. A simulated homogeneous nucleation rate of 5.0 × 10(10) cm(-3) s(-1) was obtained for our system, which is in reasonable agreement with experimental and simulation rates for homogeneous nucleation of ice at similar degrees of supercooling. This study represents our first step in a series of studies aimed at understanding the nucleation and growth of crystals of organic salts near surfaces and inside nanopores.

Reaction coordinate of incipient methane clathrate hydrate nucleation

Nucleation from solution is a ubiquitous phenomenon with relevance to myriad scientific disciplines, including pharmaceuticals, biomineralization, and disease. One prominent example is the nucleation of clathrate hydrates, multicomponent crystalline inclusion compounds relevant to the energy industry where they block pipelines and also constitute a potential vast energy resource. Despite their importance, the molecular mechanism of incipient hydrate formation remains unknown. Herein, we employ advanced molecular simulation tools (pB histogram, equilibrium path sampling) to provide a statistical-mechanical basis for extracting physical insight into the molecular steps by which clathrates form. Through testing the Mutually Coordinated Guest (MCG) order parameter, we demonstrate that both guest (methane) and host (water) structuring are crucial to accurately describe the nucleation of hydrates and determine a critical nucleus size of MCG-1 = 16 at 255 K and 500 bar. Equipped with a validated (and novel) reaction coordinate, subsequent equilibrium path sampling simulations yield the free energy barrier and nucleation rate. The resulting quantitative nucleation process is described by the MCG clustering mechanism. This constitutes a significant advance in the field of hydrates research, as the fitness of a molecular descriptor has never been statistically verified. More broadly, this work has significance to a wide range of multicomponent nucleation contexts wherein the formation mechanism depends on contributions from both solute and solvent.

Computer simulation of epitaxial nucleation of a crystal on a crystalline surface

We present results of computer simulations of crystal nucleation on a crystalline surface, in the Lennard-Jones model. Motivated by the pioneering work of Turnbull and Vonnegut [Ind. Eng. Chem. 44, 1292 (1952)], we investigate the effects of a mismatch between the surface lattice constant and that of the bulk nucleating crystal. We find that the nucleation rate is maximum close to, but not exactly at, zero mismatch. The offset is due to the finite size of the nucleus. In agreement with a number of experiments, we find that even for large mismatches of 10% or more, the formation of the crystal can be epitaxial, meaning that the crystals that nucleate have a fixed orientation with respect to the surface lattice. However, nucleation is not always epitaxial, and loss of epitaxy does affect how the rate varies with mismatch. The surface lattice strongly influences the nucleation rate. We show that the epitaxy observed in our simulations can be predicted using calculations of the potential energy between the surface and the first layer of the nucleating crystal, in the spirit of simple approaches such as that of Hillier and Ward [Phys. Rev. B 54, 14037 (1996)].

Reaction coordinates, one-dimensional Smoluchowski equations, and a test for dynamical self-consistency

We propose a method for identifying accurate reaction coordinates among a set of trial coordinates. The method applies to special cases where motion along the reaction coordinate follows a one-dimensional Smoluchowski equation. In these cases the reaction coordinate can predict its own short-time dynamical evolution, i.e., the dynamics projected from multiple dimensions onto the reaction coordinate depend only on the reaction coordinate itself. To test whether this property holds, we project an ensemble of short trajectory swarms onto trial coordinates and compare projections of individual swarms to projections of the ensemble of swarms. The comparison, quantified by the Kullback-Leibler divergence, is numerically performed for each isosurface of each trial coordinate. The ensemble of short dynamical trajectories is generated only once by sampling along an initial order parameter. The initial order parameter should separate the reactants and products with a free energy barrier, and distributions on isosurfaces of the initial parameter should be unimodal. The method is illustrated for three model free energy landscapes with anisotropic diffusion. Where exact coordinates can be obtained from Kramers-Langer-Berezhkovskii-Szabo theory, results from the new method agree with the exact results. We also examine characteristics of systems where the proposed method fails. We show how dynamical self-consistency is related (through the Chapman-Kolmogorov equation) to the earlier isocommittor criterion, which is based on longer paths.

Polymorph specific RMSD local order parameters for molecular crystals and nuclei: α-, β-, and γ-glycine

Crystal nucleation is important for many processes including pharmaceutical crystallization, biomineralization, and material synthesis. The progression of structural changes which occur during crystal nucleation are often described using order parameters. Polymorph specific order parameters have been developed for crystallization of spherically symmetric particles; however, polymorph specific order parameters for molecular crystals remain a challenge. We introduce template based polymorph specific order parameters for molecular crystals. For each molecule in a simulation, we compute the root mean squared deviation (RMSD) between the local environment around the molecule and a template of the perfect crystal structure for each polymorph. The RMSD order parameters can clearly distinguish the α-, β-, and γ-glycine polymorph crystal structures in the bulk crystal and also in solvated crystallites. Surface melting of glycine crystallites in supersaturated aqueous solution is explored using the newly developed order parameters. The solvated α-glycine crystallite has a thinner surface melted layer than the γ-glycine crystallite. α-glycine forms first out of aqueous solution, so surface melted layer thickness may provide insight into interfacial energy and polymorph selection.