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

Abstract

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.

Figure 1

Histogram of q_lm and W_l for 13 atom icosahedral, fcc, and hcp clusters, as well for 15 atom bcc and 7 atom sc cluster. Figure reproduced from ref (30) with permission. Copyright 1983 APS.

Figure 2

Distribution of number of connections per particles for liquid, bcc, and fcc structure. Figure reprinted from ref (31) with permission. Copyright 1996 AIP.

Figure 3

(a) Probability distributions of q̅₄ (solid lines) and q₄ (dashed lines) for the fcc, bcc, and hcp crystals and for the undercooled liquid (LIQ) in the Lennard-Jones system. (b) Probability distributions of q̅₆ (solid lines) and q̅₆ (dashed lines) for the same phases. Figure reprinted from ref (31) with permission. Copyright 1996 AIP.

Figure 4

Free energy surface when only Q4 and Q6 were used as CVs. Figure reprinted from ref (34) with permission. Copyright 2008 AIP.

Figure 5

Local (q6, q4, and w4) and the average (q̅6, q̅4, and w̅4) bond OPs of various solid forms are plotted for comparison. Figure reprinted from ref (35) with permission. Copyright 2014 AIP.

Figure 6

(a) Free energy landscape for the system. (b) Projection of the same landscape in the Q6-U space. Figure reprinted from ref (36) with permission. Copyright 2022 Springer.

Figure 7

(a) The particles positions with respect to a central atom (black sphere) in the reference crystal environment χ₀ and in an instantaneous configuration χ are shown in gray and orange spheres, respectively. The configuration was extracted from a trajectory of bcc sodium simulated at T = 300 K and P = 1 bar. (b) Distributions of the kernel k̅_χ0(χ) (eq 11)) for the liquid (blue) and bcc (orange) sodium at 375 K and 1 atm pressure. Figure reprinted from ref (41) with permission. Copyright 2019 AIP.

Figure 8

(a) A “point molecule” representation of glycine: r is the center of mass (COM) of the molecule, q (quaternion) indicates the absolute orientation of the molecule in the reference coordinate frame (bottom left), ψ is the (N–C–C–O) dihedral angle. (b) Construction of the pair distribution function: r is the distance vector along the COM–COM positions of the two glycine molecules projected on the first molecular coordinate frame q₁, and ψ₁ and ψ₂ are the internal degrees of freedom that are included in the pair distribution function’s definition, eq 13. Figure reprinted from ref (46) with permission. Copyright 2011 Springer.

Figure 9

(a) Molecular vectors along C=O and N–N directions in a single urea molecule. (b) A representative urea crystal structure (polymorph I) extracted from well tempered metadynamics (WTMetaD) simulation trajectory in ref (50). Figure reprinted from ref (50) with permission. Copyright 2015 Elsevier.

Figure 10

Gray and red ellipses represent molecules in the (SME) and in a crystal template, respectively. The core molecule in the template and the marked atom in the SME almost perfectly coincide after step (ii). The closest molecules in the SME are given the template molecules in step (iii). The closest of the two template molecules is preserved if two template molecules are assigned to the same SME molecule (as indicated by the dotted line). Unmatched SME molecules are also not kept if no template molecule is associated with a particular SME molecule. The unshaded template and SME molecules are “deleted” as a consequence of these criteria prior to the RMSD reduction in step (iv). Figure reprinted from ref (57) with permission. Copyright 2011 AIP.

Figure 11

MetaD simulation snapshots of low density water (LDW) structures and ice structures (I_c = cubic ice, I_sd = stacking disorder ice). On the right free-energy landscape is shown plotted using CV1 and CV2, where color pallet shows free energy in kJ/mol. Figure reprinted from ref (58) with permission. Copyright 2011 Nature.

Figure 12

FES projected on the s_H and s_S variables for (a) Na at 350 K and (b) Al at 800 K. Figure reprinted from ref (40) with permission. Copyright 2017 AIP.

Figure 13

g(r, θ) for the liquid and polymorph I of urea at 450 K with snapshots of system in both phases. Figure reprinted from ref (62) with permission. Copyright 2018 National Academy of Sciences.

Figure 14

MetaD simulation of the paracetamol system obtained biasing both KL^r̂ and Γ^rv. Figure reprinted from ref (67) with permission. Copyright 2018 ACS.

Figure 15

Example of the two unique hydrogen bonds between molecules in the form 1 crystal. Figure reprinted from ref (67) with permission. Copyright 2018 ACS.

Figure 16

Nucleation trajectory obtained from a multiple walkers metadynamics simulation biasing the positional and orientational OPs, KL_c^{r̂,r̂_OO,r̂_ON} and KL_c^v₁,v₂. Figure reprinted from ref (67) with permission. Copyright 2018 ACS.

Figure 17

Simulated XRD patterns for β-cristobalite and liquid silica at 2,400 K with system containing 1536 atoms. Figure reprinted from ref (72) with permission. Copyright 2018 National Academy of Sciences.

Figure 18

Distribution of the local structure factor S_i(q₁) in the liquid and the solid phase. Figure reprinted from ref (76) with permission. Copyright 2018 APS.

Figure 19

FES from 100 ns MetaD simulation of a 512-atom system (NaCl), using CN and V as CVs, at T = 300 K and P = 20 GPa. Figure reprinted from ref (83) with permission. Copyright 2021 APS.

Figure 20

A comparison among three CV profiles: (a) XRD intensity CV s⁰¹¹ for Na, (b) HLDA CV s^H, (c) VAC CV s^V for Na. Figure reprinted from ref (78) with permission. Copyright 2019 AIP.

Figure 21

(a) Evolution of χ₆ profile with time, (b) reweighted free energy profile of WTMetaD. Polymorph I (in yellow) and IV (in green) in (a,b) figure. Multiple transitions are visible from Polymorph I to IV. Error bar is shown in shaded blue color. Figure reprinted from ref (91) with permission. Copyright 2021 ACS.

Figure 22

Path-CV, f(Q(r)) as a function of d-AFED (top) and MetaD (bottom) simulation time. The red and blue spheres indicate Mo atoms in bcc and A15 phases, respectively. Figure reprinted from ref (93) with permission. Copyright 2019 AIP.

Figure 23

(a) S(k) plot of CO₂ crystal is shown with respective Miller indices. (b) The evolution of CVs profile with time obtained from OPES MetaD. The color is based on value of the Miller index (111) which is the first descriptor of Deep-LDA. Red color represents the crystal, and blue color represents the liquid, (c) Free energy profile of CO₂ crystallization, and (d) crystal and liquid structure of CO₂ are shown. Figure reprinted from ref (97) with permission. Copyright 2021 Taylor and Francis.

Figure 24

Comparative study between Deep-LDA and Deep-TICA driven simulation: Time evolution CV profile of (A) Deep-LDA and (B) Deep-TICA. (C) and (D) indicate the correlation between Deep-LDA and Deep-TICA CVs with the fraction of diamond atoms. White circles represent the average values of two CVs in liquid and solid phases while the dotted gray lines interpolate between them. Figure reprinted from ref (98) with permission. Copyright 2021 National Academy of Sciences.