Computer Simulations of Homogeneous Nucleation of Benzene from the Melt

Theoretical Investigation into Polymorphic Transformation between β-HMX and δ-HMX by Finite Temperature String

Polymorphic transformation is important in chemical industries, in particular, in those involving explosive molecular crystals. However, due to simulating challenges in the rare event method and collective variables, understanding the transformation mechanism of molecular crystals with a complex structure at the molecular level is poor. In this work, with the constructed order parameters (OPs) and K-means clustering algorithm, the potential of mean force (PMF) along the minimum free-energy path connecting β-HMX and δ-HMX was calculated by the finite temperature string method in the collective variables (SMCV), the free-energy profile and nucleation kinetics were obtained by Markovian milestoning with Voronoi tessellations, and the temperature effect on nucleation was also clarified. The barriers of transformation were affected by the finite-size effects. The configuration with the lower potential barrier in the PMF corresponded to the critical nucleus. The time and free-energy barrier of the polymorphic transformation were reduced as the temperature increased, which was explained by the pre-exponential factor and nucleation rate. Thus, the polymorphic transformation of HMX could be controlled by the temperatures, as is consistent with previous experimental results. Finally, the HMX polymorph dependency of the impact sensitivity was discussed. This work provides an effective way to reveal the polymorphic transformation of the molecular crystal with a cyclic molecular structure, and further to prepare the desired explosive by controlling the transformation temperature.

Machine Learning Classification of Local Environments in Molecular Crystals

Identifying local structural motifs and packing patterns of molecular solids is a challenging task for both simulation and experiment. We demonstrate two novel approaches to characterize local environments in different polymorphs of molecular crystals using learning models that employ either flexibly learned or handcrafted molecular representations. In the first case, we follow our earlier work on graph learning in molecular crystals, deploying an atomistic graph convolutional network combined with molecule-wise aggregation to enable per-molecule environmental classification. For the second model, we develop a new set of descriptors based on symmetry functions combined with a point-vector representation of the molecules, encoding information about the positions and relative orientations of the molecule. We demonstrate very high classification accuracy for both approaches on urea and nicotinamide crystal polymorphs and practical applications to the analysis of dynamical trajectory data for nanocrystals and solid-solid interfaces. Both architectures are applicable to a wide range of molecules and diverse topologies, providing an essential step in the exploration of complex condensed matter phenomena.

Finite Temperature String with Order Parameter as Collective Variables for Molecular Crystal: A Case of Polymorphic Transformation of TNT under External Electric Field

An external electric field is an effective tool to induce the polymorphic transformation of molecular crystals, which is important practically in the chemical, material, and energy storage industries. However, the understanding of this mechanism is poor at the molecular level. In this work, two types of order parameters (OPs) were constructed for the molecular crystal based on the intermolecular distance, bond orientation, and molecular orientation. Using the K-means clustering algorithm for the sampling of OPs based on the Euclidean distance and density weight, the polymorphic transformation of TNT was investigated using a finite temperature string (FTS) under external electric fields. The potential of mean force (PMF) was obtained, and the essence of the polymorphic transformation between o-TNT and m-TNT was revealed, which verified the effectiveness of the FTS method based on K-means clustering to OPs. The differences in PMFs between the o-TNT and transition state were decreased under external electric fields in comparison with those in no field. The fields parallel to the c-axis obviously affected the difference in PMF, and the relationship between the changes in PMFs and field strengths was found. Although the external electric field did not promote the convergence, the time of the polymorphic transformation was reduced under the external electric field in comparison to its absence. Moreover, under the external electric field, the polymorphic transformation from o-TNT to m-TNT occurred while that from m-TNT to o-TNT was prevented, which was explained by the dipole moment of molecule, relative permittivity, chemical potential difference, nucleation work and nucleation rate. This confirmed that the polymorphic transformation orientation of the molecular crystal could be controlled by the external electric field. This work provides an effective way to explore the polymorphic transformation of the molecular crystals at a molecular level, and it is useful to control the production process and improve the performance of energetic materials by using the external electric fields.

Experimental and Theoretical Investigation of External Electric-Field-Induced Crystallization of TKX-50 from Solution by Finite-Temperature String with Order Parameters as Collective Variables for Ionic Crystals

External electric fields are an effective tool to induce phase transformations. The crystallization of ionic crystals from solution is a common phase transformation. However, understanding of mechanisms is poor at the molecular level. In this work, we carried out an experimental and theoretical investigation of the external electric-field-induced crystallization of TKX-50 from saturated formic acid solution by finite-temperature string (FTS) with order parameters (OPs) as collective variables for ionic crystals. The minimum-free-energy path was sketched by the string method in collective variables. The results show that the K-means clustering algorithm based on Euclidean distance and density weights can be used for enhanced sampling of the OPs in external electric-field-induced crystallization of ionic crystal from solution, which improves the conventional FTS. The crystallization from solution is a process of surface-mediated nucleation. The external electric field can accelerate the evolution of the string and decrease the difference in the potential of mean forces between the crystal and the transition state. Due to the significant change in OPs induced by the external electric field in nucleation, the crystalline quality was enhanced, which explains the experimental results that the external electric field enhanced the density, detonation velocity, and detonation pressure of TKX-50. This work provides an effective way to explore the crystallization of ionic crystals from solution at the molecular level, and it is useful for improving the properties of ionic crystal explosives by using external electric fields.

Finite temperature string by K-means clustering sampling with order parameters as collective variables for molecular crystals: application to polymorphic transformation between β-CL-20 and ε-CL-20

Polymorphic transformation of molecular crystals is a fundamental phase transition process, and it is important practically in the chemical, material, biopharmaceutical, and energy storage industries. However, understanding of the transformation mechanism at the molecular level is poor due to the extreme simulating challenges in enhanced sampling and formulating order parameters (OPs) as the collective variables that can distinguish polymorphs with quite similar and complicated structures so as to describe the reaction coordinate. In this work, two kinds of OPs for CL-20 were constructed by the bond distances, bond orientations and relative orientations. A K-means clustering algorithm based on the Euclidean distance and sample weight was used to smooth the initial finite temperature string (FTS), and the minimum free energy path connecting β-CL-20 and ε-CL-20 was sketched by the string method in collective variables, and the free energy profile along the path and the nucleation kinetics were obtained by Markovian milestoning with Voronoi tessellations. In comparison with the average-based sampling, the K-means clustering algorithm provided an improved convergence rate of FTS. The simulation of transformation was independent of OP types but was affected greatly by finite-size effects. A surface-mediated local nucleation mechanism was confirmed and the configuration located at the shoulder of potential of mean force, rather than overall maximum, was confirmed to be the critical nucleus formed by the cooperative effect of the intermolecular interactions. This work provides an effective way to explore the polymorphic transformation of caged molecular crystals at the molecular level.

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.

Collective Variables for Crystallization Simulations-from Early Developments to Recent Advances

Crystallization is an important physicochemical process which has relevance in material science, biology, and the environment. Decades of experimental and theoretical efforts have been made to understand this fundamental symmetry-breaking transition. While experiments provide equilibrium structures and shapes of crystals, they are limited to unraveling how molecules aggregate to form crystal nuclei that subsequently transform into bulk crystals. Computer simulations, mainly molecular dynamics (MD), can provide such microscopic details during the early stage of a crystallization event. Crystallization is a rare event that takes place in time scales much longer than a typical equilibrium MD simulation can sample. This inadequate sampling of the MD method can be easily circumvented by the use of enhanced sampling (ES) simulations. In most of the ES methods, the fluctuations of a system's slow degrees of freedom, called collective variables (CVs), are enhanced by applying a bias potential. This transforms the system from one state to the other within a short time scale. The most crucial part of such CV-based ES methods is to find suitable CVs, which often needs intuition and several trial-and-error optimization steps. Over the years, a plethora of CVs has been developed and applied in the study of crystallization. In this review, we provide a brief overview of CVs that have been developed and used in ES simulations to study crystallization from melt or solution. These CVs can be categorized mainly into four types: (i) spherical particle-based, (ii) molecular template-based, (iii) physical property-based, and (iv) CVs obtained from dimensionality reduction techniques. We present the context-based evolution of CVs, discuss the current challenges, and propose future directions to further develop effective CVs for the study of crystallization of complex systems.

Stepwise Homogeneous Melting of Benzene Phase I at High Pressure

We investigate, using molecular dynamics simulations, the spontaneous homogeneous melting of benzene phase I under a high pressure of 1.0 GPa. We find an apparent stepwise transition via a metastable crystal phase, unlike the direct melting observed at ambient pressure. The transition to the metastable phase is achieved by rotational motions, without the diffusion of the center of mass of benzene. The metastable crystal completely occupies the whole space and maintains its structure for at least several picoseconds, so that the phase seems to have a local free energy minimum. The unit cell is found to be unique—no such crystalline structure has been reported so far. Furthermore, we discuss the influence of pressure control on the melting behavior.

Computational Study on Homogeneous Melting of Benzene Phase I

Molecular-dynamics simulations are used for examining the microscopic details of the homogeneous melting of benzene phase I. The equilibrium melting temperatures of our model were initially determined using the direct-coexistence method. Homogeneous melting at a higher temperature is achieved by heating a defect- and surfacefree crystal. The temperature-dependent potential energy and lattice parameters do not indicate a premelting phase even under superheated conditions. Further, statistical analyses using induction times computed from 200 melting trajectories were conducted, denoting that the homogeneous melting of benzene occurs stochastically, and that there is no intermediate transient state between the crystal and liquid phases. Additionally, the critical nucleus size is estimated using the seeding approach, along with the local bond order parameter. We found that the large diffusive motion arising from defect migration or neighbor-molecule swapping is of little importance during nucleation. Instead, the orientational disorder activated using the flipping motion of the benzene plane results in the melting nucleus.

Solubility curves and nucleation rates from molecular dynamics for polymorph prediction – moving beyond lattice energy minimization

Seeded nucleation simulations allow for the screening of low energy structures to predict final structure present in solution.

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.

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.

Metadynamics studies of crystal nucleation

Crystallization processes are characterized by activated events and long timescales. These characteristics prevent standard molecular dynamics techniques from being efficiently used for the direct investigation of processes such as nucleation. This short review provides an overview on the use of metadynamics, a state-of-the-art enhanced sampling technique, for the simulation of phase transitions involving the production of a crystalline solid. In particular the principles of metadynamics are outlined, several order parameters are described that have been or could be used in conjunction with metadynamics to sample nucleation events and then an overview is given of recent metadynamics results in the field of crystal nucleation.

Bézier curve string method for the study of rare events in complex chemical systems

We present a new string method for finding the most probable transition pathway and optimal reaction coordinate in complex chemical systems. Our approach evolves an analytic parametric curve, known as a Bézier curve, to the most probable transition path between metastable regions in configuration space. In addition, we demonstrate that the geometric properties of the Bézier curve can be used to construct the optimal reaction coordinate near the most probable reaction path, and can further be used to devise a ranking vector capable of identifying precisely which collective variables are most important for governing the transition between metastable states. We discuss the algorithmic details of the Bézier curve string method, analyze its stability, accuracy and efficiency, and illustrate its capabilities using model potential energy functions. In particular, we use the degree elevation property of Bézier curves to develop an algorithm that adaptively learns the degree polynomial necessary to accurately represent the most probable transition path. Subsequently, we apply our method to the isomerization of alanine dipeptide, and demonstrate that the reaction coordinate obtained from the Bézier curve string method is in excellent agreement with the optimal reaction coordinate constructed from an aimless shooting and maximum likelihood procedure. Finally, we apply our method to a large complex system and study the homogenous nucleation of benzene from the melt. In these two examples, we illustrate that the ranking vector correctly identifies which collective variables govern these chemical transitions.

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.