111${\dagger}$ Authors have contributed to the review equally.

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

Neha† Department of Chemistry, Indian Institute of Technology, Delhi,
Hauz Khas, New Delhi 110016, India  Vikas Tiwari† Department of Chemistry,
Indian Institute of Technology, Delhi,
Hauz Khas, New Delhi 110016, India  Soumya Mondal† Department of Chemistry,
Indian Institute of Technology, Delhi,
Hauz Khas, New Delhi 110016, India  Nisha Kumari† Department of Chemistry,
Indian Institute of Technology, Delhi,
Hauz Khas, New Delhi 110016, India  Tarak Karmakar  Department of Chemistry,
Indian Institute of Technology, Delhi,
Hauz Khas, New Delhi 110016, India<EMAIL_ADDRESS>
###### Abstract

Crystallization is one of the most important physicochemical processes 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 timescales 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. An ES method enhances the fluctuations of a
system’s slow degrees of freedom, called collective variables (CVs), by
applying a bias potential, and thereby transforms the system from one state to
the other within a short timescale. The most crucial part of an ES method 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.

##  1 1\. INTRODUCTION

Understanding phase transition using computer simulations has been a prime
focus of the simulators - starting from the early method developers to the the
current community. Alongside the enrichment of our fundamental understanding
of the symmetry-breaking transition from a non-symmetric liquid state to a
symmetric crystalline phase, this process is of great interest and bears
relevance to the pharmaceutical industry as well as the environment.
Investigation of the crystallization mechanism, optimizing crystallization
conditions, predicting crystal shape, and obtaining thermodynamics and
kinetics of crystal growth/dissolution processes require significant time
investment as well as monetary resources.

Computer simulations, especially the molecular dynamics (MD) methods, have
been probably the most convenient way to delineate microscopic details of the
early stage of the crystallization process and calculate the thermodynamics
and kinetics of the process. Unfortunately, in the context of computer
simulations, crystallization is a rare event that, in most cases, takes place
in a timescale ranging from milliseconds to seconds. The brute-force MD
simulations are limited by short timescales in the range of nano- to
microseconds that are inadequate to study crystallization. To circumvent this
issue, several enhanced sampling (ES) simulations methods have been developed
over the years. 1, 2, 3, 4, 5, 6, 7, 8 The central aspect of most of these
methods is to define one or more variables that are functions of atomic
coordinates and describe the system’s slow degrees of freedom. These variables
are called order parameters (OPs) or collective variables (CVs) (note that
there are subtle differences between these two nomenclatures). In ES
simulations, the fluctuations of these CVs are enhanced to sample the
metastable states. In the context of crystallization, the CVs are designed
such a way that they can distinguish between the particles in the crystal and
liquid or dispersed in solution phases.

The CVs that are routinely used in crystallization simulations can be majorly
divided into four categories, (i) spherical particle-based CVs such as the
most popular Steinhardt parameters that are useful in simulating spherical
particles (atomic, metallic, and colloids systems), (ii) molecular CVs
(template-based, root mean square deviation-based, local crystallinity) that
are used to crystallize molecular systems in pre-defined structures, (iii)
physical property-based CVs such as volume/density, pair correlation
functions, structure factor and X-ray diffraction peaks, and entropy-enthalpy,
and (iv) CVs that are derived from linear and non-linear (machine learning)
dimensionality reduction of basic OPs (distance, coordination number, angles,
etc.).

In this review, we provide a systematic overview of the CVs that are used in
ES simulations of crystallization. In doing so, we might miss many other
important developments including those CVs that are used as classifiers
(fingerprints) to characterise disordered and various crystal polymorphs.
Discussion on these topics can be found in the literature. 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20 We conclude this review with a discussion on the
current challenges and future directions in the development of efficient CVs
for the study of crystallization of complex molecular systems.

##  2 2\. COLLECTIVE VARIABLES (CVs)

###  2.1 2.1 Spherical-particle-based

####  2.1.1 2.1.1. Steinhardt parameters

In 1983, Steinhardt et al. proposed bond-orientational order parameters (OPs)
to characterize and distinguish between solid and liquid states. 21 These OPs
identify the system’s states by measuring the symmetries of the clusters
formed during simulation. In their approach, a central atom, $i$ is considered
which forms ‘virtual bonds’ with its neighbors that are found within the
radial distance of $1.2r_{0}$ around it$i$, where $r_{0}$ is the minimum
distance in the Lennard-Jones (LJ) potential. Here these ‘virtual bonds’ do
not imply any chemical bonds rather they are imaginary lines that connect the
central atom with its neighbors. The atoms’ connectivity is defined by
spherical harmonic function, and the OP, $\bar{q}_{lm}(\vec{r})$ for particle
$i$ is defined by taking an average of these spherical harmonics over a
suitable set of neighbors around the particle $i$,

$\bar{q}_{lm}(i)=\frac{1}{N_{b}(i)}\sum_{j=1}^{N_{b}(i)}Y_{lm}(\theta(r_{ij}),\phi(r_{ij}))$
(1)

Here $N_{b}(i)$ is number of virtual bonds of particles $i$ with its
neighbors, $Y_{lm}(\theta,\phi)$ are the spherical harmonics, and
$\theta(r_{ij})$ and $\phi(r_{ij})$ are polar angles made by bonds with
respect to some reference coordinate system. The value of
$Y_{lm}(\theta,\phi)$ depends upon two intergers $l$ and $m$ and for a given
value of l there are $2l+1$ values of $m$. Hence, $q_{lm}(i)$ will have
different values for a particular value of $l$. Therefore a particular
structure will have different values of $q_{lm}(i)$ for a particular value of
$m$ for a given $l$. Therefore rotationally invariant combinations of bond
order parameters are defined as:

$\displaystyle q_{l}(i)$
$\displaystyle=\sqrt{\frac{4\pi}{2l+1}\sum_{m=-l}^{l}|{\bar{q}_{lm}(i)}|^{2}}$
$\displaystyle\hat{W_{l}}$
$\displaystyle=\frac{W_{l}}{(\sum_{l=-m}^{m}|\bar{q}_{lm}|^{2})^{\frac{3}{2}}}$

where,

$W_{l}=\sum_{\begin{subarray}{c}m_{1},m_{2},m_{3}\\\
m_{1}+m_{2}+m_{3}=0\end{subarray}}\begin{pmatrix}l&l&l\\\
m_{1}&m_{2}&m_{3}\end{pmatrix}\bar{q}_{lm_{1}}\bar{q}_{lm_{2}}\bar{q}_{lm_{3}}$
(2)

The matrix $\begin{pmatrix}l&l&l\\\ m_{1}&m_{2}&m_{3}\end{pmatrix}$ is called
Wigner $3j$ symbols. The $q_{lm}$ and $W_{l}$ are called as quadratic and
third order invariants, respectively. Figure 1 shows the histogram of $q_{lm}$
and $W_{l}$ for five clusters at $l=2,4,6,8,10$. Non-zero averages appear at
$l\geq 4$ for clusters having hcp and cubic symmetry. For icosahedral
clusters, non-zero averages occur at $l=6$. All the values are calculated for
clusters corresponding to the unit cell.

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. (Reproduced from
Ref. 21, with the permission of APS Publisher.)

In 1996, ten Wolde et al. 22 studied the crystallization of LJ system using MD
simulations. They have computed nucleation barrier, nucleation rate, and
identified pre-critical, critical and post-critical nuclei for LJ system. They
have calculated the free energy surface for nucleation by using the umbrella
sampling method.1 In order to calculate the nucleation barrier, they defined
an OP (reaction coordinate) which measures the degree of crystallinity of a
system during phase transition. They found that the local OPs introduced by
Steinhardt (Eq. 1) has almost the same values in both solid and liquid phases.
Therefore, they introduced the global Steinhardt parameters, $q_{6}$ which
vanishes in liquid and have high values in the crystalline phase. The
generalized global orientational OP, $q_{lm}$ is defined as:

$\bar{q}_{lm}=\frac{\sum_{i=1}^{N}N_{b}(i)\bar{q}_{lm}(i)}{\sum_{i=1}^{N}N_{b}(i)}$
(3)

In $\bar{q}_{lm}$, the average is taken over all $N$ particles present in the
system. The values of global orientational OPs for different crystal systems
can be seen in the Table 2.1.1,

Table 1: Global orientational order parameters for fcc, hcp, bcc, sc, and
icosahedral structures. (Reprinted with the permission of AIP publisher from
ref. 22)

|  $Q_{4}$ |  $Q_{6}$ |  $\hat{W}_{4}$ |  $\hat{W}_{6}$
---|---|---|---|---
fcc |  0.191 |  0.575 |  -0.159 |  -0.013
hcp |  0.097 |  0.485 |  0.134 |  -0.012
bcc |  0.036 |  0.511 |  0.159 |  0.013
sc |  0.764 |  0.354 |  0.159 |  0.013
Icosahedral |  0 |  0.663 |  0 |  -0.170
liquid |  0 |  0 |  0 |  0

Subsequently, ten Wolde, Frenkel and coworkers 23 introduced the scalar
product of the normalized bond orientational vectors, $\bf{q_{l}(i)}$ and
$\bf{q_{l}(j)}$ between neighboring particles $i$ and $j$, used to
differentiate between solid and liquid clusters as described below :

$\bf{q_{l}(i)}.\bf{q_{l}(j)}=\sum_{m=-l}^{l}\bf{q_{l}(i)}.\bf{q_{l}(j)}^{*}$
(4)

Here, $\bf{q_{l}(j)^{*}}$ is the complex conjugate of $\bf{q_{l}(j)}$, and
$\bf{q_{l}(i)}.\bf{q_{l}(i)}$ is equal to 1. Two neighbors, $i$ and $j$ are
connected if the dot product is greater than a threshold, say 0.5. According
to this criterion all particles in solid state are found to be connected to
each other. However this is not a sufficient condition to call a cluster as
solid-like or liquid-like because the liquid-like cluster can also be
connected frequently because of the presence of local order in liquids.
Therefore an additional condition has been included that if the number of
“connections” are above some threshold, say 6 or 8 then the particles are
solid-like, and if less, then they are liquid-like. Using this criterion, the
solid-like particles can be distinguished from the liquid-like particles as
the former will have more number of connections (coordination) than the
latter. Distribution of number of connections per particle in LJ system for
liquid, bcc, and fcc structures is shown in Figure 2.

Figure 2: Distribution of number of connections per particles for liquid, bcc,
and fcc structure. (Reprinted with the permission of AIP publisher from Ref.
22)

The approach of Steinhardt, ten Wolde, Frenkel, and coworkers has been further
extended by Eslami et al. 24 who have defined a local OP as follows,

$q_{l}(i)=\frac{1}{N_{b}(i)}\sum_{j\in
N_{b}(i)}\bf{q_{l}(i)}.\bf{q_{l}(j)}^{*}$ (5)

The $q_{l}(i)$ OP includes contributions from the first and second
coordination shell neighbors of the particle $i$ averaged over all its
neighbors.

Lenchner and Dellago proposed a variant of the Steinhardt OPs which is more
accurate to determine specific crystal structures. 25. They introduced the OPs
that are obtained by averaging the bond orientation orders over all neighbors,

$\bar{q}_{l}(i)=\sqrt{\frac{4\pi}{2l+1}\sum_{m=-l}^{l}|\bar{q}_{lm}(i)|^{2}}$
(6)

where $\bar{q}_{lm}(i)$ is defined as,

$\bar{q}_{lm}(i)=\frac{1}{\tilde{N}_{b}(i)}\sum_{k=0}^{\tilde{N}_{b}(i)}q_{lm}(k)$
(7)

where the index $k$ goes from 0 to all $N_{b}(i)$ neighbors of $i$ including
itself. The advantage of this definition are - these OPs are not restricted to
including only the first coordination shell but they can take into account the
second shell neighbors, and the solidlike and liquidlike particles can be
distinguished better due to the decreased overlap between the OPs
distributions belonging to the two states.

They have calculated the average of the probability distributions of
$\bar{q}_{4}$ and $\bar{q}_{6}$ for fcc, bcc, and hcp crystals in undercooled
liquid in which particles are interacting via LJ and Gaussian potential (Table
2.1.1 and 2.1.1). Due to the averaging procedure the overlap between order
parameter distributions in different phases decreases. As a result sharper
distinction between different phases is obtained.

Table 2: Global orientational order parameters for fcc, hcp, bcc, sc, and
icosahedral structure.( Reprinted with the permission from AIP publisher from
Ref. 25)

|  $q_{4}$ |  $\bar{q}_{4}$ |  $q_{6}$ |  $\bar{q}_{6}$
---|---|---|---|---
bcc |  0.089 988 |  0.033 406 |  -0.440 526 |  0.408 018
fcc |  0.170880 |  0.158180 |  0.507298 |  0.491385
hcp |  0.107923 |  0.084052 |  0.445384 |  0.421813
liq |  0.109049 |  0.031246 |  0.360012 |  0.161962

Table 3: Average of the distributions of q4 , $\bar{q}4$, q6 , and $\bar{q}6$
for the bcc, fcc,and hcp crystals in system interacting via Gaussian
potential.(Reprinted with the permission of AIP publisher from Ref. 25)

|  $q_{4}$ |  ${\bar{q}}_{4}$ |  $q_{6}$ |  ${\bar{q}_{6}}$
---|---|---|---|---
bcc |  0.085581 |  0.031728 |  0.437129 |  0.407515
fcc |  0.155336 |  0.134388 |  0.474079 |  0.447782
hcp |  0.109723 |  0.073369 |  0.424627 |  0.385720
liq |  0.126950 |  0.040297 |  0.375121 |  0.158913

In 2014, Tang et al. have used as CVs the local and average Steinhardt OPs
with or without supercell parameters included in the CVs definition to predict
the polymorphism of Xenon crystal at high temperature and pressure.26 As
expected, fcc and bcc structures were obtained when $Q4$ and $Q6$ were used as
CVs (see Fig 3), and along with it, the new structures including fcc with hcp
stacking faults were also obtained when the supercell parameters was included
in the CVs definition.

Figure 3: Free energy surface when only $Q4$ and $Q6$ were used as
CVs.(Reprinted with the permission of AIP publisher from Ref. 26)

The structures obtained by using local and average order parameters as CV were
compared (see Fig. 4). The fcc structure with stacking faults show two types
of structures with the local OPs were used as CV, and the fcc structure with
stacking faults were found to split into more types of structures when the
average OP was used as CV.

Figure 4: The local (q6, q4, and w4) and the average ($\bar{q}6$, $\bar{q}4$,
and $\bar{w}4$) bond OPs of various solid forms are plotted for
comparison.Reprinted with the permission of AIP publisher from Ref. 26

Very recently, Rozanov et al. has studied the phase transition of a LJ system
from its metastable liquid to crystalline state using MetaD simulations. 27
The $Q6$ along with the systems potential energy ($U$) were used as CVs. The
free energy landscape associated with this phase transformation process is
shown in Figure 5.

Figure 5:  (a)Free energy landscape for the system (b) Projection of the same
landscape in the Q6-U space.(Reprinted with the permission of Springer
publisher from Ref. 27)

The free energy barrier of crystallization at constant pressure obtained from
the MetaD simulations is in good agreement with those obtained from
experiments. This manifests the CVs effectiveness in sampling the actual
nucleation process and reproducing the experimental observations. 27

###  2.2 2.2. SOAP Kernel

####  2.2.1 2.2.1. SOAP Kernel method

A crucial property of a variable representing atomic environments is its
invariance with respect to the basic symmetries like rotation, reflection,
translation, and permutation of atoms. Steinhardt OPs are one of the
invariants used to describe atomic environments. But these OPs have certain
limitations. 28 The qualitative trend of the Steinhardt OPs is influenced by
the choice of neighborhood, and due to the discrete nature of neighborhood
definition, the neighborhood of a particle is not a continuous function of
particles coordinates. This discontinuous nature of $q_{l}$ leads to the lack
of robustness of these OPs as structure metrics. Bartok et al. has shown that
the descriptors like Steinhardt parameters are the special cases of some
general approach in which the atomic environment is defined by neighbourhood
density. 29 This approach is called Smooth Overlap Of Atomic Positions (SOAP).
In the SOAP Kernel method, each atom in a given environment is defined as the
sum of Gaussian functions centered on the neighborhood of an atom and
including that atom itself. These parameters fulfil the criteria of being
invariant and continuous functions of atomic coordinates. The SOAP kernel and
its variants have been used in many applications for identifying crystal
structures and polymorphs. 30, 31

####  2.2.2 2.2.2 Environment Similarity CV

Piaggi and Parrinello utilized a reduced definition of the SOAP kernel
approach to design a CV, called environment similarity CV (Env-CV). 32 This
method is based on the local ordering of $n$ neighbors ($i$) around a central
atom ($\bf{r}$) in a crystalline environment (Fig. 6). As done in the SOAP
kernel approach, the local density $\rho_{\chi}(\bf r)$ of the central atom in
a given environment $\chi$ is written as the sum of Gaussian functions,

$\begin{split}\rho_{\chi}(\bf r)&=\sum_{i\in\chi}e^{-|{\bf r_{i}}-{\bf
r}|^{2}/2\sigma^{2}}\end{split}$ (8)

where ${\bf r_{i}}$’s are the coordinates of the neighbors relative to the
central atom, and $\sigma^{2}$ is the variance of the Gaussian functions. The
$n$ nearest neighbors positions $\\{\bf r_{j}^{0}\\}$ of the central atom in a
reference crystal environment ($\chi_{0}$) are chosen. The difference between
the two environments $\chi$ and $\chi_{0}$ is obtained from the following
integral,

$\begin{split}k_{\chi_{0}}(\chi)&=\int dr\rho_{\chi}(\bf r)\rho_{\chi_{0}}(\bf
r)\end{split}$ (9)

where $\rho_{\chi_{0}}(\bf r)$ is the local density of the atom in the
reference crystal environment ($\chi_{0}$).

Unlike in the SOAP’s actual definition 29, in this approach, only the spatial
part has been considered. This results in a simple analytical expression of
the CV in the form of a kernel function, $k_{\chi_{0}}(\chi)$ which can be
calculated efficiently.

$\begin{split}k_{\chi_{0}}(\chi)&=\sum_{i\in\chi}\sum_{j\in\chi_{0}}\pi^{3/2}\sigma^{3}e^{-|{\bf
r_{i}}-{\bf r^{0}_{j}}|^{2}/4\sigma^{2}}\end{split}$ (10)

The kernel function in Eq. (10) is then normalized such that similarities
between identical environments $\bar{k}_{\chi_{0}}(\chi_{0})$ is equal to one.

$\begin{split}\bar{k}{{}_{\chi_{0}}(\chi)}&=\frac{k_{\chi}(\chi_{0})}{k_{\chi_{0}}(\chi_{0}))}\\\
&=\frac{1}{n}\sum_{i\in\chi}\sum_{j\in\chi_{0}}e^{-|{\bf r_{i}}-{\bf
r^{0}_{j}}|^{2}/4\sigma^{2}}\end{split}$ (11)

In a system of $N$ particles, for each particle ($i$ = 1,..,$N$) the kernel
function, $\bar{k}_{\chi_{0}}(\chi_{i})$ is calculated using Eq. (11). The
particles having $\bar{k}_{\chi_{0}}(\chi_{i})$ $>$ ${k_{0}}$, where
${k_{0}}$=0.5, are classified as crystalline or liquid by a continuous and
differentiable switching function ($s^{O}_{i}$) as follows,

$\begin{split}s_{i}&=1-\frac{1-(\bar{k}_{\chi_{0}}(\chi_{i})/k_{0})^{p}}{1-(\bar{k}_{\chi_{0}}(\chi_{i})/k_{0})^{q}}\end{split}$
(12)

The variable $s_{i}$ has values in the range from 0 to 1; for atoms in the
solution $s_{i}$ $\approx$ 0 while those in a perfect crystalline environment,
$s_{i}$ $\approx$ 1 (Fig. 6(b)). The parameters $p$ and $q$ control the
steepness and the range of the switching function. The kernel function defined
in Eq.(10) is in the similar spirit as that of the local metric order by
Martelli, Car, and co-workers. 33

Figure 6: (a) The particles positions with respect to a central atom (black
sphere) in the reference crystal environment, $\chi_{0}$ and in an
instantaneous configuration $\chi$ are shown in grey 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
$\bar{k}{{}_{\chi_{0}}(\chi)}$ (eq. 11) ) for the liquid (blue) and bcc
(orange) sodium at 375 K and 1 atm pressure.(Reprinted with the permission of
AIP publisher from Ref. 32)

The Env-CV has been successfully used in the study of phase diagram of sodium
and aluminium using a multithermal-multibaric enhanced sampling simulation
approach. 32 The non-rotationally invariant nature of the CV facilitates the
crystallization of a defect-free crystal. Niu et al. extended the application
of this CV to calculate the phase diagram of Gallium modelled using a DeepNN
potential. 34 Recently, Piaggi et al used the Env-CV to study molecular
system, in particular, ice nucleation from water modelled using a deep
learning based potential model. The Env-CV’s use was not restricted to single-
component systems crystallizations. Karmakar et al. used this CV to nucleate
NaCl from its supersaturated aqueous solution using a combination of MetaD and
constant chemical potential MD approaches. 35

###  2.3 2.3. Molecular ordering

####  2.3.1 2.3.1. A generalized set of OPs

For highly symmetric systems consisting atoms or spherical particles
(colloids), Steinhardt bond OPs are found effective. However, for complex low
symmetric systems extending such OPs are way difficult to execute.

To define OPs for complex crystals, preliminary information can be generated
from a normal unbiased MD simulation of the crystal. An OP to define complex
crystallization is to take advantage of the structural properties of the given
crystal or consider the atomic coordinates relevant to a specific molecular
arrangement in a crystal. A new method to design an OP to study complex
crystal systems is presented here.36

To build OPs for molecular crystals, a generalized pair correlation function
consisting of all relevant variables that represent the crystal structure is
introduced by Santiso et al. 37 Before moving ahead with pair correlation
functions, it is important to understand the idea of point molecule
representation (Fig. 7). In point molecule representation, the crystal system
is reduced to defining - (i) position (the center of mass of each molecule in
the crystal) (ii) absolute orientation (containing a set of molecule-centered
coordinate axes), and (iii) the internal degrees of freedom (define the
internal structure of the original molecule). The OPs are then defined as the
product of the probability density functions $f$ as:

$G(r,q,{\psi},{\psi}^{{}^{\prime}})\approx\sum_{\alpha=1}^{\infty}f_{\alpha}^{r}(r)f_{\alpha}^{\hat{r}}{(\hat{r})}f_{\alpha}^{q}(q)f_{\alpha}^{\psi}(\psi,\psi^{{}^{\prime}})$
(13)

Figure 7: (a) A ‘point molecule’ representation of glycine - $r$ is the centre
of mass (COM) of the molecule, $q$ (quaternion) indicates the absolute
orientation of the molecule in the reference coordinate frame (bottom left),
$\psi$ is the (N-C-C-O) dihedral angle. (b) Construction of the pair
distribution function - $\bf{r}$ distance vector along the COM-COM positions
of the two glycine molecules projected on the first molecular coordinate frame
$q1$, and $\psi_{1}$ and $\psi_{2}$ are the internal degrees of freedom that
are included in the pair distribution function’s definition, Eq.13. (Reprinted
with the permission of Springer publisher from Ref. 37)

where $r$ is the distance between the center of mass of two molecules,
${\hat{r}}$ is the bond orientation, and the relative orientation of first
molecule with respect to the second molecule is represented by $q$, and
${\psi}^{{}^{\prime}}$ is the relative internal degrees of freedom with
respect to ${\psi}$ (the internal degrees of freedom for the first molecule).
Now, we can define the pair distribution function as the probability that a
molecule has an internal degrees of freedom between ${\psi}$ and ${\psi}$ \+
$d{\psi}$ with a neighbor with internal degrees of freedom between
${\psi}^{{}^{\prime}}$ and $\psi^{{}^{\prime}}$ \+ $d{\psi}^{{}^{\prime}}$,
relative orientation between $q$ and $q+dq$ at a position between $r$ and
$r+dr$ with respect to the first molecule. The values $r$, $q$, ${\psi}$ and
${\psi}^{{}^{\prime}}$ uniquely represent the peaks in the pair distribution
function and define the crystal structure. In order to define an OP, there is
a need to choose models for each function appearing in the probability density
function, $G(r,q,\psi,\psi^{{}^{\prime}})$ in Eq. 13 (Fig.7(b)). Parameters
are estimated for the models using an unbiased simulation. The distribution of
distance between center of mass of the molecules, $f_{\alpha}^{r}(r)$ is
approximated using a Gaussian function,

$f_{\alpha}^{r}(r)\approx\frac{1}{\sqrt{2\pi}\sigma_{\alpha}}exp[-\frac{(r-r_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}]$
(14)

Where $\sigma_{\alpha}$ is the standard deviation, $\alpha$ is the peak
corresponding to the mean center of mass distance r$\alpha$. For each of the
peaks in the pair distribution function, the distribution of bond orientation,
$f_{\alpha}^{(\hat{r})}(\hat{r})$ is approximated using the Fisher
distribution,

$f_{\alpha}^{(\hat{r})}(\hat{r})\approx\frac{\kappa_{\alpha}}{sinh\kappa_{\alpha}}exp(\kappa_{\alpha}\hat{r}_{\alpha}^{T}\hat{r})$
(15)

Where $\kappa_{\alpha}$ is the concentration parameter, and $\alpha$ is the
peak corresponding to the mean bond orientation $\hat{r}\alpha$. For the
relative orientation distribution, $f_{\alpha}^{q}(q)$, 4D Bingham
distribution can be used, as the orientation vectors are directionless on a 4D
unit sphere. However, whatever studies have been done so far, bipolar Watson
distribution can be a good approximation around each peaks in the pair
distribution function $G(r,q,\psi,\psi^{,})$:

$f_{\alpha}^{q}(q)\approx\frac{1}{{}_{1}F_{1}(1/2,2,\xi_{\alpha})}exp[\xi_{\alpha}(q_{\alpha}.q)^{2}]$
(16)

where $\xi_{\alpha}$ is the concentration parameter, $q\alpha$ is the mean
relative orientation, $1F1(1/2,2,x)$ is the confluent hypergeometric function,
and the 4D dot product is denoted by the ‘.’ symbol. Finally, the accurate
model for the internal degrees of freedom is chosen on a case by case basis.
When describing a crystal structure, the internal degrees of freedom that are
considered are atom distances, angles, and dihedrals. For distance between two
atoms, a Gaussian distribution model is used, and for angles and dihedrals,
the von Mises distribution is generally used.

$\phi_{\alpha,i}^{d}=\sum_{j\neq
i}\frac{1}{\sqrt{2\pi}\sigma_{\alpha}}exp[-\frac{(r-r_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}]$
(17)

where $i$ is the the central molecule, and $r_{ij}$ is the center-of-mass
separation between molecule $i$ and its neighbors, $j$. Similar to this, one
can define OPs that consider both bond distances and bond orientations,

$\phi_{\alpha,i}^{bo}=\frac{1}{\sqrt{2\pi}\sigma_{\alpha}}\frac{\kappa_{\alpha}}{sinh\kappa_{\alpha}}\sum_{j\neq
i}exp[-\frac{(r-r_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}]exp[\kappa_{\alpha}\hat{r}_{\alpha}^{T}\hat{r}]$
(18)

where $\hat{r}_{ij}$ represents the bond orientation vector projected onto the
frame with the molecule $i$ at its centre. Additional order criteria that are
sensitive to relative orientations of molecules include,

$\phi_{\alpha,i}^{ro}=\frac{1}{\sqrt{2\pi}\sigma_{\alpha}}\frac{1}{{}_{1}F_{1}(1/2,2,\xi_{\alpha})}\times\sum_{j\neq
i}exp[-\frac{(r-r_{\alpha})^{2}}{2\sigma_{\alpha}^{2}}]exp[\xi_{\alpha}(q_{\alpha}.q)^{2}]$
(19)

where, $F_{1}(1/2,2,\xi_{\alpha})$ is the confluent hypergeometric function,
$q_{\alpha}$ is the relative orientation, $\xi_{\alpha}$ is a concentration
parameter, and ‘.’ indicates a 4D dot product. Likewise, OPs that take into
account a molecule’s internal configuration can be defined.

The “local” or per-molecule and per-peak, OPs mentioned above can be used to
quantify which molecule initiates the process of ordering as well as its
extent. However, it is impractical for use in complex system, instead, a
“global” OP obtained by adding up either or both of the indices $i$ and
$\alpha$ in Eqs. (17)-(19) would be more convenient to use. The design and
application of this order parameters is presented in the study of (i)
crystallization of $\alpha$-glycine from solution which demonstrates how to
cope with a nonsymmetric molecule with adaptable internal degrees of freedom,
(ii) the crystallization (nucleation) 38 of benzene from the melt, which
serves as an example of how the OPs for a relatively high symmetric molecule
are constructed, and (iii) solid-solid polymorph transformation of
terephthalic acid.

####  2.3.2 2.3.2 Local crystallinity order

Giberti, Salvalaglio, and Parrinello have developed CVs based on local
crystallinity order for studying molecular crystals. 39, 40, 41, 42, 43 The
synergistic effect of local density fluctuations and molecular orientational
ordering has been embedded in the CVs definition. The advantages of this CV
are two-fold: (i) It allows the system to explore the free energy surface
without the prior knowledge about the global crystal symmetry, and (ii) its
use can be extended to crystallization in multicomponent systems such as
solution crystallization. 40, 43

The global CV is defined as the sum of individual molecules crystallinity
order ($\Gamma_{i}$) in a system containing $N$ molecules,

$s=\frac{1}{N}\sum_{i=1}^{N}\Gamma_{i}$ (20)

Local density of a central $i^{th}$ molecule is calculated by calculating its
coordination number (CN) with respect to its neighbors, $j$. From the central
atom, distances ($r_{ij}$) from its $j$ neighbors within a cut-off distance,
$r_{cut}$ are calculated. To decide weather $j^{th}$ particle will be
considered as neighbor of $i^{th}$ particle a switching function (Fermi),
$f(r_{ij})$ is used,

$f(r_{ij})=\frac{1}{1+e^{a(r_{ij}-r_{cut})}}$

The summation of $f(r_{ij})$ defines the neighbor density as:

$n_{i}=\sum_{j}f(r_{ij})$ (21)

Another switching function, $\rho_{i}$ is defined to calculate the local
density as a function of coordination number as:

$\rho_{i}=\frac{1}{1+e^{-}b(n_{i}-n_{cut})}$ (22)

where $a$ and $b$ are used to tune the slope of exponential functions in the
switching functions.

Figure 8: (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 WTMetaD simulation trajectory in ref. 41.(Reprinted from Ref. 41,
Copyright (2015), with permission from Elsevier.)

To calculate the orientation between molecules $i$ and $j$, a function
$\theta_{ij}$ which is a function of the angle between two molecular vectors
is defined (Fig. 8(a)). Due to the effect of temperature the orientation
between molecules fluctuates around an average value ($\bar{\theta}_{k}$).
Therefore the fluctuation is expressed in terms of Gaussian functions centered
around the ($\bar{\theta}_{k}$) between two molecular vectors,

$\theta_{ij}=\sum_{k=1}^{k_{max}}e^{-\frac{(\theta_{ij}-\bar{\theta_{k}})^{2}}{2\sigma_{k}^{2}}}$
(23)

Here $k_{max}$ is the maximum number of angles that define the local molecular
orientations. Based on these spatial and local molecular orientational OPs,
the molecular crystallinity is defined as,

$\Gamma_{i}=\frac{\rho_{i}}{n_{i}}\sum_{j=1}^{N}f_{ij}\Theta_{ij}$ (24)

A CV based on the total number of crystalline molecule is obtained by summing
over the individual crystallinity orders,

$s=\sum_{i=1}^{N}\Gamma_{i}$ (25)

In ref. 41, the fraction of molecules that are in crystalline environment has
been used as a global CV,

$s=\frac{1}{N}\sum_{i=1}^{N}\left[\frac{\rho_{i}}{n_{i}}\sum_{j=1}^{N}f_{ij}\sum_{k=1}^{k_{max}}e^{-\frac{(\theta_{ij}-\bar{\theta_{k}})^{2}}{2\sigma_{k}^{2}}}\right]$
(26)

The molecular crystallinity CV has been used in my recent studies either to
characterize crystal-like molecules in a multi-component system or to
construct the bias potential in WTMetaD simulations studying crystallization
from melt or solution. 40, 43, 44, 45

####  2.3.3 2.3.3. Molecular RMSD

Inspired by the shape matching OP by Keys et al. 46 and the pattern matching
variables by Shetty et al. 47, Duff and Peters introduced the template based
polymorph specific OP which computes root mean square deviation (RMSD) between
a tagged crystal molecule in a simulation system and a molecular template in a
perfect crystal.48 The RMSD-based OP can audibly differentiate between
different polymorphs of a crystal structure.

Though theoretically comparable to the OPs proposed by Shetty et al., this
order parameter is computationally simpler. By using a template of the crystal
structure and the environment of a tagged molecule in a simulation, the method
calculates the desired RMSDs. For the purpose of showcasing the new technique,
Peters et al. carried out the same order parameter diagnostics as Lechner and
Dellago25 here for molecular crystal polymorphs. They showed that their
technique is capable of differentiating between the bulk crystal structures of
the three glycine polymorphs without any overlap in the order parameter
distributions. Additionally, in solvated glycine crystallites, the $\alpha$,
$\beta$, and $\gamma$ glycine polymorph structures may be distinguished using
the local RMSD based OPs. The approach offers a broad framework that makes it
simple to design OPs for a range of molecular crystals.

The local molecular order is obtained by matching a tagged solute in the
crystal and its simulated microenvironment (SME) comprise of neighboring
solutes to its corresponding “central” template molecule and its neighbors
present in a perfect crystal (from Cambridge Crystallographic Database or a
modelled equilibrium lattice structure). If the crystal has $m$ solutes, there
will be $m$ distinct templates to build the crystal structure. If one deals
with $k$ number of crystal polymorphs, for each polymorph ($k$) one needs to
define $m_{k}$ molecular templates. Fig. 9 shows the template-matching
procedure.

Figure 9: The grey 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).(Reprinted with the permission of AIP publisher from
Ref. 48)

Here we briefly mention the steps involved in the RMSD OPs development. To
reduce the computational cost, the hydrogen atoms are not considered in the
RMSD calculations, and the functional groups that can adopt degenerate
configurations are replaced with non-degenerate surrogate functional groups.
After this ‘molecular pruning’ step, the RMSD OP calculation is initiated. The
overall process has the following steps - (i) at first, a solute molecule is
tagged and the solutes surrounding it within a cutoff radius are used to
define the SME. (ii) This is followed by the rotation and translation of the
central molecule in the crystal template to minimize the RMSD between the
tagged and the central molecule. The same amount of rotation and translation
is performed with the entire template. (iii) Subsequently, each molecule in
the SME is matched with its nearest template molecule. (iv) Finally, all
molecules in the SME are matched with the molecules in the crystal (polymorph)
template. The final step’s reduced RMSD value serves as an OP for the tagged
molecule that is particular to unit-cell-member and polymorph. For each member
of the polymorph $k$ unit cell, repeat steps (i-iv). Then choose the $k^{th}$
polymorph’s smallest minimized RMSD value.

Both for bulk crystallites and a small crystal in solution, the RMSD OPs can
differentiate between various polymorphs with clarity. The new OPs should make
it possible to simulate small molecules, something that was previously only
achievable for supercooled simple liquids.

###  2.4 2.3. Physical-property-based

Recently, the use of CVs related to a system’s physical properties that can be
calculated experimentally has become a preferred choice for studying
crystallization using simulations. This is because, in general, the value of a
particular physical property is known experimentally for all the states of a
system, and it can be directly used to construct a useful CV. An important
aspect of using the physical property as a CV is that it does not require
prior knowledge of a system’s crystalline state. The physical properties that
have been utilized to construct a CV are radial distribution functions, XRD
peak intensities, entropy and enthalpy, and system’s volume (density).

####  2.4.1 2.3.1. Radial distribution function

Nada in 2020, proposed a method in which the radial distribution functions
(RDFs) can be utilized as a CV for the formation of water polymorphs. 49 They
performed MetaD simulations using two CVs defined as two discrete oxygen-
oxygen RDFs represented by Gaussian window functions. Different polymorphs of
ice such as cubic, stacking disorder 50 (consists of cubic and hexagonal),
high pressure ice VII, layered ice with an ice VII, and layered ice with an
unknown structure were identified from the MetaD simulation trajectory (Fig.
10).

Figure 10: 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
colour pallet shows free energy in kJ/mol.(Reprinted with the permission of
Nature from Ref. 49)

####  2.4.2 2.3.2. Entropy and Enthalpy

The CVs discussed in sections 2.1 and 2.2 were constructed based on known
crystal structures, and thus they are not effective in discovering other
possible polymorphic phases of the crystal. Hence, there was a need to
construct CVs that can sample the states without any prior knowledge of the
crystal structure. Keeping this in mind, in 2017, Piaggi et. al. proposed the
use of enthalpy and entropy surrogates as CVs. 31 This choice was based on two
simple facts - (i) ‘enthalpy and entropy’ - that do not predict any feature of
the crystal structure a priori, and (ii) there is a trade-off between
‘enthalpy’ and ‘entropy’ during the crystallization which, in turn, describes
the transitions between metastable states. Although ‘enthalpy’ is easy to
estimate, the ‘entropy’ calculation is a non-trivial task. However, in the
context of crystallization, we do not require an exact definition of entropy
to bias the system; an approximate equation involving only two body
correlations, derived from an expression where excess entropy per atom is
expressed as an infinite series of terms involving multiparticle correlation
functions suffices the need. 51

The two CVs constructed using enthalpy and entropy are defined below:

$s_{H}=\frac{U(R)+PV}{N}$ (27) $s_{S}=-2\pi\rho
k_{B}\int_{0}^{\infty}[g_{m}(r)ln(g_{m}(r))-g_{m}(r)+1]r^{2}dr$ (28)

where, $g_{m}(r)$ is the mollified version of the radial distribution function
to ensure the function’s continuity.

$g_{m}(r)=\frac{1}{4\pi N\rho r^{2}}\sum_{i\neq
j}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-(r-r_{ij})^{2}}/{2\sigma^{2}}$ (29)

where, $\sigma$ is the broadening parameter, $r_{ij}$ is the distance between
$i^{th}$ and $j^{th}$ particle, and $\rho$ is the system’s density. These CVs
were used to study the crystallization of Na and Al from their molten states
(Fig.11).

Figure 11: FES projected on the $s_{H}$ and $s_{S}$ variables for - (a) Na at
350 K and (b) Al at 800 K.(Reprinted with the permission of APS Physics from
Ref. 31)

Mendels et al. in 2018, extended the applicability of these CVs to study a
multicomponent system, silver iodide (AgI). 52 In this work, they have
predicted the existence of an $\alpha$ phase of AgI which is stabilized by
strong entropic contributions in comparison to the enthalpically-favored
$\beta$ phase.

Although, these CVs were successful in predicting polymorphism in atomic
crystals like Na and Al but they cannot be used for molecular crystals because
molecules do not have a spherical symmetry hence they can have different
orientations in space, and depending on its orientation they can exist in
different polymorphic forms. The above defined CVs do not take into account
these orientations and thus are less efficient in the case of molecular
crystals.

To tackle this problem, in 2018 Piaggi et al. proposed the use of
orientational entropy as a CV for predicting polymorphisms in molecular
crystals. 53 In this case, along with the spatial distances they included
molecular orientations ($s_{\theta}$).

$s_{\theta}=-2\pi\rho
k_{B}\int_{0}^{\infty}\int_{0}^{\pi}[g_{m}(r,\theta)ln(g_{m}(r,\theta))-g_{m}(r,\theta)+1]r^{2}sin\theta
drd\theta$ (30)

where, $\theta$ is the angle between two vectors $v_{i}$ and $v_{j}$
describing the orientation of molecules i and j.

$\theta=arccos\left(\frac{v_{i}.v_{j}}{|v_{i}||v_{j}|}\right)$ (31)

In principle, at least three angles are required to define relative
orientation of molecule in space, for example, Euler angles
$\phi,\theta,\psi$. Hence, our function would look like
$g(r,\phi,\theta,\psi)$ which is not very convenient to work with. As this
function has four variables making our simulations complicated and less
efficient. So, instead of taking one CV with 4 variables better alternative is
to take 2 CVs $s_{\theta_{1}}$ and $s_{\theta_{2}}$ defining two different
relative orientation of molecules using two angles ${\theta_{1}}$ and
${\theta_{2}}$. To understand the behaviour of $g(r,\theta)$ for liquid and
solid phases we can see an example of Urea at 450 K (Fig. 12). Here, $\theta$
represents the direction of dipole moment in urea. From the figure 12 we can
observe that the liquid exhibits some structure at very short distances
whereas, in polymorph I, a well-defined structure exists at long distances.
One of the main characteristics of polymorph I revealed by $g(r,\theta)$ is
that molecules have parallel and antiparallel dipole moments.

Figure 12: g(r, $\theta$) for the liquid and polymorph I of urea at 450 K with
snapshots of system in both phases.(Reprinted with the permission of national
academy of sciences from Ref. 53)

Further, using these CVs WTMetaD simulations were performed for urea and
naphthalene at 450 K and 300 K, respectively. These temperatures are close to
the melting temperatures of both substances. A large number of transitions to
different crystal forms have been observed. To identify and classify the
polymorphs formed during the simulation, a similarity finding strategy has
been adopted in Ref. 54 and 55. between two given configurations. The distance
between two g(r, $\theta$)’s i.e., the divergence is calculated as

$D(g_{1}\parallel
g_{2})=\int_{0}^{\infty}\int_{0}^{\pi}\left[g_{1}(r,\theta)ln\frac{g_{1}(r,\theta)}{g_{2}(r,\theta)}-g_{1}(r,\theta)+g_{2}(r,\theta)\right]\times
r^{2}sin\theta drd\theta$ (32)

This is Kullback-Leibler divergence for non-normalized functions with a minima
at $g_{1}=g_{2}$.56 As $D(g_{1}\parallel g_{2})$ is not symmetric, it cannot
be used as distance, hence distance is given by

$d(g_{1},g_{2})=\frac{D(g_{1}\parallel g_{2})+D(g_{2}\parallel g_{1})}{2}$
(33)

Using hierarchical clustering and average distance between points in two
clusters trajectory of urea was analyzed and different crystalline forms and
liquid form were successfully distinguished.

Amodea et al. used this entropy surrogate CV along with the potential energy
CV to study the effect of cooling rate during Ni3Al nanoparticle freezing. 57
In that work, they found that by adjusting the cooling rate of Ni3Al
nanoparticles one can stabilize an out of equilibrium polymorph, BCC DO$-3$
structure.

####  2.4.3 2.3.3. Information Entropy

In another work, Gobbo et al. used a CV based on the relative information
entropy along with the Santiso and Trout’s pair-distribution function based
CVs to study crystallization of benzene and paracetamol. 58 The point molecule
representation 37 which is characterized by the position of its molecular
center ($r$) and two orientation vector ($v_{1}$ and $v_{2}$) have been used.
The per-molecule OP is written as,

$\begin{split}\Gamma_{i}^{rv}=\frac{1}{n_{i}}\sum_{j\neq
i}s(|r_{i}-r_{j}|)\times\sum_{\alpha=1}^{M}\\{exp[-((|r_{i}-r_{j}|-d_{\alpha})^{2}/2\sigma_{d_{\alpha}}^{2})]\\\
\times
exp[-((\theta(v_{1i},v_{1j})-\theta_{1\alpha})^{2}/2\sigma_{\theta_{1\alpha}}^{2})]\\\
\times
exp[-((\theta(v_{2i},v_{2j})-\theta_{2\alpha})^{2}/2\sigma_{\theta_{2\alpha}}^{2})]\\}\end{split}$
(34)

where, $M$ is the number of peaks in the joint distribution of distances and
angles, $\theta_{\alpha}$ and $d_{\alpha}$ are the peak centers, and $\sigma$
is the width of the Gaussian. $s(r)$ is a switching function.

The global average of these OPs is given as

$\Gamma^{rv_{1}v_{2}}=\frac{1}{N}\sum_{i}\Gamma_{i}^{rv_{1}v_{2}}$ (35)

These CVs discussed above have some major limitations such as - they are (i)
not able to distinguish between different polymorphs, (ii) prone to have
degeneracies (different configurations giving same value of CVs), and (iii)
not able to describe more complex structures such as paracetamol. Hence to
overcome these limitations, the authors of ref. 58 used a similar approach of
utilizing entropy based CVs as first proposed in ref. 53. Here instead of
using the entropy surrogate as a CV, they used its distributions to
differentiate the ordered state from the disordered state taking into account
the long range correlations which can give better resolution of the states. To
construct the CV, relevant quantities can be selected and the relative
probability density $p$ is build on-the-fly and compared with the suitable
reference distribution, $q$. The relative entropy is calculated using the
Kullback-Leibler divergance (KLD) method59,

$KL(q\parallel p)=\int q(x)ln\left(\frac{q(x)}{p(x)}\right)dx$ (36)

From the above equation it is clear that the value of KLD can only be
positive, and it is zero only when $q(x)=p(x)$. The probability density
($p(x)$) is given by

$p(x)\approx\frac{\sum_{i}w(x_{i})g(|x-x_{i}|)}{\sum_{j}w(x_{j})}$ (37)

where the sum runs over all elements, $g$ is the normalized Gaussian function,
and $w$ is the weights. To avoid numerical instability due to very small
values of $p$, the Eq. (36) is modified as follows,

$KL\left(q\parallel\left(\frac{1}{2}p+\frac{1}{2}q\right)\right)$ (38)

Kernel density estimate (KDE)60 must be evaluated on a grid to compute
integrals numerically. Hence the OP takes the final form as,

${\sum_{m}}q(m)ln\left(\frac{q(m)}{(p(m)+q(m))/2}\right)dS$ (39)

$m$ runs over all the grid points, and $dS$ is the measure of the volume
element associated with every grid point.

Using this approach, one can construct CVs of increasing complexity that can
help understand the crystallization of molecular systems. In ref. 58, the
authors have constructed a CV, $KL^{\hat{r}}$ where a set of distance vectors
between centers of molecules is considered. For every normalized distance
vector, an $azimuthal$ angle ($\theta$) and a dihedral angle ($\phi$) are
calculated. The value of the CV, $KL^{\hat{r}}$ is calculated using Eq. (39)
using uniform distribution as the reference distribution, $q$.

Figure 13: MetaD simulation of the paracetamol system obtained biasing both
$KL^{\hat{r}}$ and $\Gamma^{rv}$.(Reprinted with the permission of ACS from
Ref. 58)

For benzene, a good separation between the liquid and crystal states was
observed, and additionally, another ordered structure $C_{2}$ (possibly a
polymorph) was observed as well. MetaD simulations gave two different pathways
of form I crystal formation when the $KL^{\hat{r}}$ CV was used alone or used
together with the pair-function based CV ($\Gamma^{rv}$). In the first
pathway, orientational ordering is followed by form I crystal formation,
whereas, in the second pathway, positional ordering is followed by transition
to form I crystal.

In the case of paracetamol, when both CVs, $KL^{\hat{r}}$ and $\Gamma^{rv}$
were biased, the system efficiently sampled multiple ordered and disordered
states (Fig. 13). A few of the ordered states resemble the form I crystal of
paracetamol however, the obtained structures showed defects and unmatched
lattice parameters.

Figure 14: Example of the two unique hydrogen bonds between molecules in the
form 1 crystal.(Reprinted with the permission of ACS from Ref. 58)

So to further improve the results, authors constructed OPs using only KLD
framework because use of pair-function OPs was not able to define complexity
of paracetamol molecules. To compensate orientational contribution which was
previously described by pair-function OPs, a new KLD based OPs were
constructed using orientational vectors of molecules namely $KL_{c}^{v_{1}}$
and $KL_{c}^{v_{2}}$. Instead of using them separately which will increase
computational cost due to multidimentionality, they were combined to give one
OP given as $KL_{c}^{v_{1},v_{2}}$ = ($KL_{c}^{v_{1}}$ \+ $KL_{c}^{v_{2}}$)/2.
To incorporated effect of hydrogen bonding during form I crystal formation
$KL_{c}^{\hat{r}}$ was modified as
$KL_{c}^{\hat{r},\hat{r}_{OO},\hat{r}_{ON}}$ = ($KL_{c}^{\hat{r}}$ \+
$KL_{c}^{\hat{r}_{OO}}$ \+ $KL_{c}^{\hat{r}_{ON}}$)/3. Where last two term
take care of hydrogen bonds formed in from I crystal according to the fig 14.
In metadynamics simulation it was observed nucleation event starts around 20
ns with the formation of dimers, after this system rapidly orders to form full
crystal in 60ns simulation (Fig. 15. After analysis it was observed that
$KL_{c}^{\hat{r}_{OO}}$ and $KL_{c}^{\hat{r}_{ON}}$) were the slowest evolving
OPs which tells that hydrogen bond formation is the rate determining step in
this nucleation process. Once the hydrogen bonds are formed in proper
direction and orientation, it stabilizes the complex leading to drive the
nucleation mechanism forward. Despite success of these CVs, no hydrogen bonds
were found to be present in the dimeric species formed at first in the
nucleation mechanism and left this mechanism unclear which leaves room for
further research to understand the complexity of nucleation mechanism.

Figure 15: Nucleation trajectory obtained from a multiple walkers metadynamics
simulation biasing the positional and orientational OPs,
$KL_{c}^{\hat{r},\hat{r}_{OO},\hat{r}_{ON}}$ and
$KL_{c}^{v_{1},v_{2}}$.(Reprinted with the permission of ACS from Ref. 58)

In 2020, Song et al. used the concept of Shannon information entropy based CVs
to predict polymorphism in 1:1 cocrystal of resorcinol and urea using
adiabatic free energy dynamics (AFED). 61

####  2.4.4 2.3.4. Structure factor and XRD

Recently, Invernizzi and Niu used the concepts of structure factor and XRD-
peaks, respectively, to design suitable CVs for the enhanced sampling
simulation of crystallization. 62, 63 One of the most important properties of
a crystal is its X-ray diffraction pattern which is easily obtained from
experiments. In an XRD experiment, the scattering intensity is derived as a
function of scattering vectors as follows

$I(\vec{Q})=\sum_{i=1}^{N}\sum_{j=1}^{N}f_{i}(Q)f_{j}(Q)e^{-i\vec{Q}.(\vec{R_{i}}-\vec{R_{j})}}$
(40)

where, $R_{i}$ and $R_{j}$ are the positions of $i^{th}$ and $j^{th}$
particle, $f(Q)$ is a function of magnitude of scattering vector (Q) known as
the scattering form factor.

In Ref. 63, the spherically averaged Debye scattering function has been used
as a CV,

$I(\vec{Q})=\sum_{i=1}^{N}\sum_{j=1}^{N}f_{i}(Q)f_{j}(Q)\frac{sin(Q.R_{ij})}{Q.R_{ij}}W(R_{ij})$
(41)

where $R_{ij}$ is the distance between the atoms $i$ and $j$. A window
function $W(R_{ij})$64-65 is used to define a soft cutoff $R_{c}$ for $R_{ij}$
to avoid numerical instability while dealing with a finite-size system.

Figure 16: Simulated XRD patterns for $\beta$-cristobalite and liquid silica
at 2,400 K with system containing 1536 atoms.(Reprinted with the permission of
ACS from Ref. 64)

In general, it is obvious to choose low-theta high intensity peaks as CVs as
they provide long-range crystalline order. Ramakrishnan-Yussouff theory of
crystallization also suggests to use highest peak of the structure factor as
freezing order parameter.66 Figure 16 shows XRD patterns for
$\beta$-cristobalite silica which shows most intense peak at {111} and {022}
and the liquid silica where the peaks loose sharp features. The two CVs,
$s_{1}$ and $s_{2}$ were therefore defined as follows,

$s_{1}=I(Q_{\\{111\\}})$ (42) $s_{2}=I(Q_{\\{022\\}})$ (43)

Using these CVs, the FES for silica crystallization was obtained. Successful
implementation of the XRD peak as CV has opened a new and better class of CVs
which can be used for studying crystallization without any prior knowledge of
the crystalline structure. Further to this development, Bonati et al. used the
local structure factor as a CV to study the nucleation of silicon from its
melt using a deep neural network potential for Si.67

The Debye formula has been modified into individual atomic contributions,

$S(q)=1+\frac{1}{N}\sum_{i=1}^{N}\sum_{j\neq
1}^{N}\frac{sin(qr_{ij})}{qr_{ij}}=\frac{1}{N}\sum_{i}^{N}S_{i}(q)$ (44)

Here, every atom is assigned to its own structure factor $S_{i}(q)$ (Fig. 17)
defined as

$S_{i}(q)=1+\sum_{j\neq 1}^{N_{n}}\frac{sin(qr_{ij})}{qr_{ij}}$ (45)

Where, the sum is over all $N_{n}$ neighbors of atom $i$ which are in a cutoff
distance of $r_{c}$.

In another work by Niu et al. in 2019, both XRD peak intensities and surrogate
of translational entropy were used to study ice nucleation from water and
their temperature dependence.68

Figure 17: Distribution of the local structure factor $S_{i}(q_{1})$ in the
liquid and the solid phase.(Reprinted with the permission of APS physics from
Ref. 67)

In their work, the CV based on scattering peak intensities was constructed
using a linear combination of seven descriptors

$s_{X}=s_{100}+s_{002}+s_{101}+\alpha(s_{100}^{xy}+s_{\bar{1}20}^{xy})+\beta
s_{002}^{xz}+\gamma s_{002}^{yz}$ (46)

here, first three peaks have high intensity, next two correspond to
intensities of two main peaks of one single honeycomb bilayer which is
projected into the XY plane, and the last two are the first main peak of the
layers that are vertical to the honeycomb bilayer in a x-z and y-z plane
respectively. $\alpha$, $\beta$ and $\gamma$ are the weights for corresponding
descriptors which has the values 2, 1 and 1 respectively in this work.

The use of XRD peaks as CVs has gained popularity in recent times for the
investigation of crystallization processes. There are many articles published
recently where XRD peaks have been utilized as CVs. 69, 70, 71, 72, 73 In
2021, Ahlawat et al. used XRD peaks as CVs to study phase transitions in
methylammonium lead iodide (MAPbI3) and formamidinium lead iodide (FAPbI3).71
They also found a low temperature crystallization pathway for the
$\alpha$-FAPbI3. In another work by Deng et al. $1^{st}$ and $2^{nd}$ XRD peak
intensities were chosen as CVs to study crystallization of silica using
enhanced sampling method and further combined with machine learning method to
find out relationship of structure and mechanical properties of silica.70 In
2021, Lodesani et al. also utilized these CVs to study crystallization path of
lithium disilicate through metadynamics simulations where they modified the
equation 41 to take only silicon atoms to calculate XRD peak intensities for
the purpose of reducing computational cost.72 In another recent work, they
used XRD peak intensity based CVs to study thermodynamics of silica
crystallization into $\beta$-cristobalite.73

####  2.4.5 2.3.5. Coordination number and Volume as CVs

More recently, Badin et al. studied pressure induced B1-B2 phase transition in
NaCl using metadynamics where they used coordination number ($CN$) and volume
($V$) as 2D collective variables. 74 The choice of $CN$ as CV was motivated by
the generic rule of high pressure chemistry which states that pressure induced
transitions are accompanied by an increase of $CN$ in the 1st coordination
sphere. Also there is a significant change in the volume of the system during
the pressure induced structural transition. In pressure induced phase
transitions, it was a smart choice to utilize very basic properties i.e., $CN$
and $V$ as the CVs to study the B1-B2 transitions in NaCl which is accompanied
by transfer of ions from the second to the first coordination shell. The
average coordination number between $Na^{+}$ and $Cl^{-}$ is calculated using
following switching function

$CN=\frac{2}{N}\sum_{\begin{subarray}{c}i\in Na^{+}\\\ j\in
Cl^{-}\end{subarray}}\left(1+\left(\frac{r_{ij}-d_{0}}{r_{0}}\right)^{6}\right)^{-1}$
(47)

where $r_{ij}$ is the distance between $i^{th}$ cation and $j^{th}$ anion and
N is the total number of of ions. $d_{0}$ and $r_{0}$ are the parameters of
the switching function which can be chosen according to the need. Using these
CVs free energy surface obtained for the B1-B2 transition of NaCl crystal
(Fig. 18).

Figure 18: 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.(Reprinted with the
permission of APS physics from Ref. 74)

###  2.5 2.4. Dimensionality reduction based

In the previous sections, we presented a large number of CVs that have been
developed and used to carry out crystallization simulations. While they are
effective on their own merits, for practical use, in enhanced sampling
simulations, only a small number of such CVs can be used. However, as Russo
and Tanaka75 pointed out, crystallization involves the ordering of multiple
OPs, and to study such a process, one needs to deal with a large number of
CVs/OPs. To alleviate this problem, various dimensionality reduction
techniques have been used that condense a large number of CVs into one or two-
dimensional ones. Here we briefly discuss some of those dimensionality methods
that have been used to design CVs for crystallization simulations.

####  2.5.1 2.4.1. Harmonic Linear Discriminant Analysis (HLDA)

Mendels et al. developed a method, Harmonic Linear Discriminant Analysis
(HLDA) to find CVs from a set of descriptors ${\bf d(R)}$ collected from
metastable states of a system. In general to construct HLDA CVs, a set of
descriptors, ${\mathbf{d}(R)}$ are calculated from unbiased simulations of a
system’s metastable phases. Subsequently, the averages and variances are used
to define the ‘between class’
(${\mathbf{S}_{b}}=({\bf\mu_{A}}-{\bf\mu_{B}})({\bf\mu_{A}}-{\bf\mu_{B}})^{T}$)
and ‘within class’ ($S_{w}=(\frac{1}{\sum_{A}}+\frac{1}{\sum_{B}})$) matrices,
respectively. The highest separation between the two states is obtained by
maximizing the Fischer’s ratio, $\frac{W^{T}S_{b}W}{W^{T}S_{w}W}$, with
respect to an $N_{d}$-dimensional projection vector, $W$. The value of $W^{*}$
that maximizes the Fischer’s ratio is obtained as,

$W^{*}=(\frac{1}{\sum_{A}}+\frac{1}{\sum_{B}})(\mu_{A}-\mu_{B})$ (48)

Finally, the HLDA CV is obtained as,

$s_{HLDA}(R)=W^{*T}{\bf d(R)}$ (49)

This method has been applied in the study chemical reaction 76 and folding of
a mini-protein. 77 Recently, Zhang et al. used this approach to find out
suitable CVs for crystallization of Na and Al from their molten states. 69
They have used a set of high intense XRD peaks (see Eq.41) of crystalline Na
($d_{1}=\tilde{I^{011}}$,$d_{2}=\tilde{I^{002}}$,$d_{3}=\tilde{I^{112}}$, and
$d_{4}=\tilde{I^{022}}$) and Al
($d_{1}=\tilde{I^{111}}$,$d_{2}=\tilde{I^{002}}$,$d_{3}=\tilde{I^{022}}$, and
$d_{4}=\tilde{I^{113}}$) as descriptors to derive the HLDA CVs. Two sets of
WTMetaD simulations were carried out - in the first, a single peak of the XRD
was biased, and in the second set, the HLDA CV, $s^{H}$. From Fig. 19, it is
clear that the HLDA CV, $s^{H}$ outperforms the single peak-based CV in
sampling the solid and liquid states.

Figure 19: A comparison among three CV profiles - (a) XRD intensity CV
$s^{011}$ for Na, (b) HLDA CV $s^{H}$, (c) VAC CV $s^{V}$ for Na. (Reprinted
with the permission of AIP from Ref. 69)

####  2.5.2 2.4.2. Time-lagged Independent Component Analysis (TICA) and
Variational Approach to Conformational Dynamics (VAC)

The time-lagged independent component analysis (TICA) linearly combines a set
of input descriptors, $\bf{d}_{k}(\bf{R_{t}})$, k = 1…$N_{d}$ to construct a
CV as, $s_{i}({\bf{R}})=\sum_{k=1}^{N_{d}}b_{ik}d_{k}$. The TICA variant
developed by Pande and Noè provides a way to optimally choose the expansion
coefficients, $\bf{b_{i}}$ by solving the eigenvalue problem,

$\bf{\tilde{C}}(\tau).\bf{b_{i}}=\bf{\tilde{C}}(0)\lambda_{i}\bf{b_{i}}$ (50)

where, $\bf{\tilde{C(0)}}$ is the covariance matrix at time 0, and
$\bf{\tilde{C(\tau)}}$ is the time lagged covariance matrix obtained as,
$C_{mn}(\tau)=\langle r_{m}(0),r_{n}(\tau)\rangle$ where,
$r_{k}(\tau)=d_{k}(\tau)-d_{k}$. $\lambda_{i}$ is the $i^{th}$ eigenvalue. The
eigenvalues are arranged in descending order, and the eigenvector having the
largest eigenvalue corresponding to the slowest degree of freedom is used as a
CV.

Usually, TICA components are obtained from a long unbiased simulation in which
the system visits metastable states multiple times as done in Refs. 78, 79,
80. McCarty and Parrinello 81 showed that a WTMetaD (biased) trajectory in
which frequent transitions between the system’s metastable states are obtained
can also be used to obtain the TICA components. However, in the latter case,
one has to obtain the scaled time from the biased simulation time,

$d\tau=e^{(-\beta V(s(\bf{R_{t}})-c(t))}dt$ (51)

Zhang et al 69 used this idea to develop a VAC CV, $s^{V}$ which was based on
the linear combination of a selected set of XRD peaks,
$I^{011},I^{002},I^{112}$, and $I^{022}$. Both HLDA CV, $s^{H}$ and the VAC
CV, $s^{V}$ exhibited improved efficiency compared to that of the single XRD
peak-based CV, $s^{011}$ (Fig. 19).

####  2.5.3 2.4.3. Spectral gap optimization of order parameters (SGOOP)

Tiwary and coworkers used the spectral gap optimization of order parameters
(SGOOP) approach to construct an one-dimensional reaction coordinate (RC) to
study nucleation of urea crystals. 82 In SGOOP 83, at first, a short MetaD
simulation is performed by taking a trial CV
($f=c_{1}\psi_{1}+c_{2}\psi_{2}+...+c_{d}\psi_{d}$) to estimate the stationary
density. Postprocessing optimization is performed in the space of mixing
coefficients $(c_{1},c_{2},..,c_{d})$ to find out best CVs with maximum
spectral gap.

Figure 20: (a) Evolution of $\chi_{6}$ profile with time, (b) Reweighted free
energy profile of WTMetaD. Polymorph I (in yellow) and IV (in green) in (a,b)
figure. Multiple transition are visible from Polymorph I to IV. Error bar is
shown in shaded blue color. (Reprinted with the permission of ACS from Ref.
82)

In ref. 82, to study urea nucleation, the RC has been defined as a linear
combination of six different OPs, entropy ($S$), enthalpy($H$), coordination
number, averaged angles $(\overline{\theta_{1}},\overline{\theta_{2}})$, and
pair orientaional entropy $(S_{\theta_{1}},S_{\theta_{2}})$.

$\chi_{6}=0.885S_{\theta_{1}}+0.330S_{\theta_{{2}}}+0.906\overline{\theta_{1}}+1.0\overline{\theta_{2}}-0.011N-0.017H$
(52)

It is clear from the coefficients of the RC defined in Eq.52, the pair
orientaional entropy $(S_{\theta_{1}},S_{\theta_{2}})$, specifically
$S_{\theta_{1}}$ has maximum weight to the RC indicating its dominant role in
nucleation events. The CV profile obtained from the WTMetaD simulations with
the SGOOP 1d-RC shows multiple transitions of the system to various metastable
states (Fig. 20(a)). These states correspond to different polymorphs of urea.
The calculated FE profile (Fig. 20(b)) indicates greater stability of
polymorph I than the polymorph-IV which is in agreement with the experiments.

####  2.5.4 2.4.4. Neural-Network-based Path Collective Variable (NN-PCV)

Rogal et al. 84 have used Behler-Parrinello symmetry functions 85, 86 and
Steinhardt parameters ($Q_{l},l=6,7,8$, see section 2.1.1) as input
descriptors for a feed-forward NN to obtain per-atom CVs ($q^{\alpha}_{i}$)
corresponding to a particular crystal structure ($\alpha=$ A15, fcc, bcc, hcp,
disordered structure). Subsequently, these atomic CVs were used to define the
global CVs as, $Q^{bcc}=(1/N)\sum_{i=1}^{N}q_{i}^{\alpha}$. Finally, the path
CV between two states, say $A15$ and $bcc$ is defined as,

$f({\bf{Q(r)}})=\frac{1}{P-1}\frac{\sum_{k=1}^{P}(P-1)exp[-\lambda\lvert{\bf{Q(r)}}-{\bf{Q_{k}(r)}}\lvert^{2}]}{\sum_{k=1}^{P}exp[-\lambda\lvert{\bf{Q(r)}}-{\bf{Q_{k}(r)}}\lvert^{2}]}$
(53)

where, $\bf{Q(r)}$ is a point on the two-dimensional $\\{Q^{A15},Q^{bcc}\\}$
space, ${\bf{Q_{k}}},k=1...P$ are the nodal points along the path of
transformation from $A15$ to $bcc$ states,
$\lvert{\bf{Q(r)}}-{\bf{Q_{k}(r)}}\lvert^{2}$ is the square distance, and
$\lambda$ is a parameter. Starting from A15 phase, the value of the Path-CV
increases from 0 to 1 reaching the bcc phase. Another Path-CV,
$z({\bf{Q(r)}})$ perpendicular to $f({\bf{Q(r)}})$, is defined to measure the
distance from the path of transformation. One can use $z({\bf{Q(r)}})$ either
to construct bias or as a restraint potential.

Figure 21: Path-CV, $f({\bf{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.(Reprinted with the permission of AIP physics
from Ref. 84)

This Path-CV was then used in d-AFED/TAMD and MetaD simulations to enhance the
phase transition between the $A15$ to $bcc$ phases of Molybdenum and calculate
the associated free energy profile (not shown).

####  2.5.5 2.4.5. Deep-LDA

So far we have discussed a few methods that are used to linearly combine a set
of descriptors to construct low dimensional CVs. Recently, a few non-linear
methods based on deep learning have been developed and used in the context of
crystallization. Deep-LDA developed Bonati et al. is one such method in which
a deep NN is appended with a LDA layer in the penultimate step of the
dimension reduction setup. 87

In the LDA method, one uses $N_{d}$ number of descriptors, d(R) which are the
functions of atom coordinates and calculate ‘within class’,
$\bf{S_{w}}=\frac{1}{2}(\sum_{S}+\sum_{L})$, and ‘between class,
$\bf{S_{b}}=(\bf{(d_{S}-d_{L})(d_{S}-d_{L})}^{T}$ matrices. Here in the
context of crystallization, the subscript $S$ and $L$ refer to the ‘Solid’ and
‘Liquid’ states, respectively. The maximum separation between the two states
is obtained by maximising the Fisher’s ratio,
$f(\bf{w})=\frac{\bf{wS_{b}w}^{T}}{\bf{wS_{w}w}^{T}}$. The value of the
projection vector, $\bf{w^{*}}$ that maximizes $f(\bf{w})$ is obtained from
the generalized eigenvalue problem, $\bf{S_{b}w_{i}}=\nu_{i}S_{w}w_{i}$. The
LDA CV finally takes the form, $s=\bf{w^{*}}^{T}d(\bf{R})$

In the Deep-LDA method, the d(R) are fed to a feed-forward deep NN which
results in an $N_{h}$ dimension output $\bf{h}$ (hidden layer). The ‘within
class’ and ‘between class’ matrices of dimension $N_{h}\times N_{h}$ are then
calculated in the $\bf{h}$ basis. The eigenvalue of the lowest eigenvector
from the Fisher’s generalized equation is used as the loss function to
optimize the NN weights. Finally, the deep-LDA CV is obtained as,
$s=\bf{w}^{T}\bf{h}$

The efficiency of a Deep-LDA CV depends on the quality of the input
descriptors. As common descriptors, one can use atom coordinates, distances,
and coordination numbers. However, most of these order parameters are
localized and include short-range orders. To study crystallization which
involves long-range ordering of molecules in a periodic lattice, Karmakar et
al. 88 used the square root of the three dimensional structure factor peaks as
input descriptors for Deep-LDA.

$d({\bf{k}})=\sqrt{s({\bf{k}})}=\frac{1}{\sqrt{N}}\bigg{\rvert}\sum_{i=1}^{N}e^{-i{\bf{k.R_{i}}}}\bigg{\rvert}$
(54)

where, $\bf{k}$ is the 3D scattering vector, and ${\bf{R_{i}}},i=1...N$ are
$N$ atoms coordinates. An appropriate choice of the $\bf{k}$ vectors results
in the formation of a particular crystal lattice. The CV by its construction
is not rotationally invariant, however, this feature favors the growth of the
crystal aligned with the MD box.

This approach has been applied in study of NaCl and CO2 crystallization from
their molten/fluid phase. In particular, for the case of CO2, the $s{\bf{k}}$
peaks with Miller indices, $(111),(012),(121),(302),(132)$, and $(004)$ were
used to describe the input descriptors (Eq. 54). The Deep-LDA CV was used in
the OPES simulations to study the phase transitions.

Figure 22: (a) S(k) plot of $CO_{2}$ 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
descriptors of Deep-LDA. Red color represent the crystal and blue color
represent the liquid, (c) Free energy profile of $CO_{2}$ crystallisation, and
(d) Crystal and liquid structure of $CO_{2}$ are shown.(Reprinted with the
permission of Taylor and Francis from Ref. 88)

A large number reversible transitions between the solid to liquid phase is
clearly visible from Fig. 22(b) manifesting the effectiveness of the Deep-LDA
CV. Similar sampling efficiency has been observed in the case of NaCl
crystallization (Fig. not shown). In both cases, the Deep-LDA-based CVs gave
better separation among the different states than in the single peak-based
CVs. Deep-LDA based CVs give new route to efficient study of crystallization
and find out the stability of crystal phase relative to its liquid phase.

####  2.5.6 2.4.6. DeepTICA

In section 2.4.2, we have discussed TICA and VAC approaches to construct CVs
from unbiased and biased simulation trajectories having multiple transitions
between metastable states. Recently, Bonati et al. extended the VAC approach
and developed its non-linear variant using a deep NN. 89

The NN of Deep-TICA takes a set of descriptors, $\bf{d(R_{t}})$ and
$\bf{d(R_{t+\tau}})$ as input features and returns a set of latent variables,
$\mathbf{h}_{\theta}(\bf{d(R_{t})}$ and
$\mathbf{h}_{\theta}(\bf{d(R_{t+\tau})}$. Subsequently, the covariance
matrices are calculated using these latent variables,

$\displaystyle C_{mn}(0)$ $\displaystyle=\langle h_{m}(t)h_{n}(t)\rangle$
$\displaystyle C_{mn}(\tau)$ $\displaystyle=\langle
h_{m}(t+\tau)h_{n}(t+\tau)\rangle$

The eigenvalues ($\lambda_{i}$) are then obtained from the generalized
eigenvalue equation,

$\bf{\tilde{C}}(\tau).\bf{b_{i}}=\bf{\tilde{C}}(0)\lambda_{i}\bf{b_{i}}$ (55)

The Deep-NN is optimized by minimizing the loss function ($\mathcal{L})$
defined as a sum of the first $N$ eigenvalues,

$\mathcal{L}=-\sum_{i=1}^{N}\lambda_{i}(\theta)^{2}$ (56)

In this way, the Deep-TICA network provides as output a set of eigenfunctions
that are used as Deep-TICA CVs in OPES simulations.

Figure 23: 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 grey lines interpolate between them.
(Reprinted with the permission of national academy of sciences from Ref. 89)

This approach has been applied to the study of prototypical processes -
alanine dipeptide conformational dynamics, folding of a mini protein
(chignolin), and crystallization of liquid Si. To fit the context of this
review, here we discuss only the case of Si crystallization. A Deep-LDA CV was
developed from a set of three-dimensional structure factor peaks, $S(\bf{k})$
(discussed in Section 2.3.4), and used in OPES simulations to sample the phase
transitions. From this trajectory, the Deep-TICA CVs were obtained. Compared
to the Deep-LDA CV, the Deep-TICA CV exhibited improved sampling between the
two states (Fig. 23) manifesting the importance of incorporation of dynamical
information in the development of an efficient CV.

##  3 DISCUSSION

In this review, we discussed some of the important order parameters or
collective variables that have been developed and used in enhanced sampling
simulations to study crystallization. The early OPs developments were based on
spherical particles that were mostly used to study atomic or metallic systems.
Attempts have been made to extend their application in complex multi-component
materials and molecular systems. However, later studies revealed that the
spherical particle-based OPs - the Steinhardt’s parameters and their variants
are not sufficient to capture molecular orientation in periodic crystals. The
atoms’ density fluctuations alone cannot fully describe the crystallinity
order. The development of local molecular OPs provided a leap toward this
goal, and these CVs have been used to sample the nucleation and growth of
organic crystals. So far, most of the systems in which these CVs were tested
consist of small, mostly rigid organic molecules. Their application in large
flexible organic crystals is, however, scarce. This is due to the fact that
for a large flexible molecular system, one needs to define a large number of
CVs that describe molecular ordering. In ES simulation, one cannot use so many
CVs, and in fact, the use of more than three variables is already a tedious
task. Dimensionality reduction-based methods are useful in such a context. In
this review, we have briefly discussed a few linear and non-linear (NN)
methods that can help construct low-dimensional CVs from a large set of
descriptors. It is important to note that the success of any dimensionality
reduction-based CV development method relies on the quality of the input
descriptors, and we have seen in a few examples that the linear combinations
of either XRD peak intensities 69, 90, 71 or entropy-based descriptors 82 are
found effective, while in another set of examples, the three-dimensional
structure factor-based CVs were non-linearly combined to develop efficient NN
CVs. 88, 89

Despite the enormous success of the above-mentioned approaches, the
development of effective CVs for large fluxional molecules such as active
pharmaceutical ingredients and biomolecules (peptides) is far from reality.
The large conformational space intrinsic to these systems makes it challenging
to design efficient CVs. ML-based dimensionality reduction methods 84, 87, 88,
89 with RMSDs as effective descriptors can be tested for such systems. Among
other possibilities, one can combine local atoms contacts combined with
systems properties such as configurational entropy and enthalpy as possible
CVs. The multicanonical approaches 32, 8, 89 along with NN-based CVs 91, 84,
85, 92, 19, 87, 89 may open up new possibilities for the study of phase
transitions of complex systems.

6\. AUTHOR INFORMATION
Corresponding Author*
Tarak Karmakar \- Department of Chemistry, Indian Institute of Technology,
Delhi, Hauz Khas, New Delhi - 110016; ORCID: 0000-0002-8721-6247; E-mail:
<EMAIL_ADDRESS>

Authors
Neha \- Department of Chemistry, Indian Institute of Technology, Delhi, Hauz
Khas, New Delhi - 110016
Vikas Tiwari \- Department of Chemistry, Indian Institute of Technology,
Delhi, Hauz Khas, New Delhi - 110016
Soumya Mondal \- Department of Chemistry, Indian Institute of Technology,
Delhi, Hauz Khas, New Delhi - 110016
Nisha Kumari \- Department of Chemistry, Indian Institute of Technology,
Delhi, Hauz Khas, New Delhi - 110016

Notes:
The authors declare no competing financial interest.
† Neha, VT, SM, and NK have contributed equally to the review article.

7\. ACKNOWLEDGEMENTS
Neha, SM, NK thank IIT Delhi for institute Ph. D. fellowship, and VT
acknowledges Prime Minister Research Fellowship. TK thanks IIT Delhi seed
grant for financial support.

## References

  * Torrie and Valleau 1997 Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. _Journal of Computational Physics_ 1997, _23_ , 187–199
  * Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping free-energy minima. _Proc. Natl. Acad. Sci. U. S. A._ 2002, _99_ , 12562–12566
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. _Phys. Rev. Lett._ 2008, _100_ , 020603
  * Valsson and Parrinello 2014 Valsson, O.; Parrinello, M. Variational approach to enhanced sampling and free energy calculations. _Phys. Rev. Lett._ 2014, _113_ , 090601
  * Valsson et al. 2016 Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint. _Annu. Rev. Phys. Chem._ 2016, _67_ , 159–184
  * Invernizzi and Parrinello 2020 Invernizzi, M.; Parrinello, M. Rethinking Metadynamics: from bias potentials to probability distributions. _J Phys Chem Lett_ 2020, _11_ , 2731–2736
  * Debnath and Parrinello 2020 Debnath, J.; Parrinello, M. Gaussian mixture-based enhanced sampling for statics and dynamics. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 5076–5080
  * Invernizzi et al. 2020 Invernizzi, M.; Piaggi, P. M.; Parrinello, M. Unified Approach to Enhanced Sampling. _Phys. Rev. X_ 2020, _10_ , 041034
  * Chau and Hardwick 1998 Chau, P.-L.; Hardwick, A. A new order parameter for tetrahedral configurations. _Molecular Physics_ 1998, _93_ , 511–518
  * Errington et al. 2003 Errington, J. R.; Debenedetti, P. G.; Torquato, S. Quantification of order in the Lennard-Jones system. _The Journal of chemical physics_ 2003, _118_ , 2256–2263
  * Radhakrishnan and Trout 2002 Radhakrishnan, R.; Trout, B. L. A new approach for studying nucleation phenomena using molecular simulations: Application to CO 2 hydrate clathrates. _The Journal of chemical physics_ 2002, _117_ , 1786–1796
  * Hawtin et al. 2008 Hawtin, R. W.; Quigley, D.; Rodger, P. M. Gas hydrate nucleation and cage formation at a water/methane interface. _Physical Chemistry Chemical Physics_ 2008, _10_ , 4853–4864
  * Peters 2009 Peters, B. Competing nucleation pathways in a mixture of oppositely charged colloids: Out-of-equilibrium nucleation revisited. _The Journal of chemical physics_ 2009, _131_ , 244103
  * Angioletti-Uberti et al. 2010 Angioletti-Uberti, S.; Ceriotti, M.; Lee, P. D.; Finnis, M. W. Solid-liquid interface free energy through metadynamics simulations. _Physical Review B_ 2010, _81_ , 125416
  * Geiger and Dellago 2013 Geiger, P.; Dellago, C. Neural networks for local structure detection in polymorphic systems. _The Journal of chemical physics_ 2013, _139_ , 164105
  * Long and Ferguson 2014 Long, A. W.; Ferguson, A. L. Nonlinear machine learning of patchy colloid self-assembly pathways and mechanisms. _The Journal of Physical Chemistry B_ 2014, _118_ , 4228–4244
  * Cheng et al. 2015 Cheng, B.; Tribello, G. A.; Ceriotti, M. Solid-liquid interfacial free energy out of equilibrium. _Physical Review B_ 2015, _92_ , 180102
  * Lee et al. 2019 Lee, S.; Teich, E. G.; Engel, M.; Glotzer, S. C. Entropic colloidal crystallization pathways via fluid–fluid transitions and multidimensional prenucleation motifs. _Proceedings of the National Academy of Sciences_ 2019, _116_ , 14843–14851
  * Fulford et al. 2019 Fulford, M.; Salvalaglio, M.; Molteni, C. DeepIce: A deep neural network approach to identify ice and water molecules. _Journal of Chemical Information and Modeling_ 2019, _59_ , 2141–2149
  * DeFever et al. 2019 DeFever, R. S.; Targonski, C.; Hall, S. W.; Smith, M. C.; Sarupria, S. A generalized deep learning approach for local structure identification in molecular simulations. _Chemical science_ 2019, _10_ , 7503–7515
  * Steinhardt et al. 1983 Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bond-orientational order in liquids and glasses. _Physical Review B_ 1983, _28_ , 784
  * Rein ten Wolde et al. 1996 Rein ten Wolde, P.; Ruiz-Montero, M. J.; Frenkel, D. Numerical calculation of the rate of crystal nucleation in a Lennard-Jones system at moderate undercooling. _J. Chem. Phys._ 1996, _104_ , 9932–9947
  * Ten Wolde et al. 1995 Ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. Numerical evidence for bcc ordering at the surface of a critical fcc nucleus. _Physical review letters_ 1995, _75_ , 2714
  * Eslami et al. 2017 Eslami, H.; Khanjari, N.; Muller-Plathe, F. A local order parameter-based method for simulation of free energy barriers in crystal nucleation. _Journal of Chemical Theory and Computation_ 2017, _13_ , 1307–1316
  * Lechner and Dellago 2008 Lechner, W.; Dellago, C. Accurate determination of crystal structures based on averaged local bond order parameters. _The Journal of chemical physics_ 2008, _129_ , 114707
  * Yu et al. 2014 Yu, T.; Chen, P. Y.; Chen, M.; Samanta, A.; Vanden-Eijnden, E.; Tuckerman, M. Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions. _The Journal of chemical physics_ 2014, _140_ , 06B603
  * Rozanov et al. 2022 Rozanov, E.; Protsenko, S.; Baidakov, V. Study of the Activation Barrier of Crystallization of a Metastable Liquid Using Metadynamics. _Physics of the Solid State_ 2022, _64_ , 22–25
  * Mickel et al. 2013 Mickel, W.; Kapfer, S. C.; Schröder-Turk, G. E.; Mecke, K. Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter. _Journal of Chemical Physics_ 2013, _138_ , 044501
  * Bartók et al. 2013 Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. _Physical Review B_ 2013, _87_ , 184115
  * De et al. 2016 De, S.; Bartók, A. P.; Csányi, G.; Ceriotti, M. Comparing molecules and solids across structural and alchemical space. _Physical Chemistry Chemical Physics_ 2016, _18_ , 13754–13769
  * Piaggi et al. 2017 Piaggi, P. M.; Valsson, O.; Parrinello, M. Enhancing Entropy and Enthalpy Fluctuations to Drive Crystallization in Atomistic Simulations. _Phys. Rev. Lett._ 2017, _119_ , 15701
  * Piaggi and Parrinello 2019 Piaggi, P. M.; Parrinello, M. Calculation of phase diagrams in the multithermal-multibaric ensemble. _The Journal of chemical physics_ 2019, _150_ , 244119
  * Martelli et al. 2018 Martelli, F.; Ko, H.-Y.; Oğuz, E. C.; Car, R. Local-order metric for condensed-phase environments. _Physical Review B_ 2018, _97_ , 064105
  * Niu et al. 2020 Niu, H.; Bonati, L.; Piaggi, P. M.; Parrinello, M. Ab initio phase diagram and nucleation of gallium. _Nat. Commun._ 2020, _11_ , 1–9
  * Karmakar et al. 2019 Karmakar, T.; Piaggi, P. M.; Parrinello, M. Molecular dynamics simulations of crystal nucleation from solution at constant chemical potential. _Journal of chemical theory and computation_ 2019, _15_ , 6923–6930
  * Murugan 2005 Murugan, N. A. Orientational Melting and Reorientational Motion in a Cubane Molecular Crystal: A Molecular Simulation Study. _The Journal of Physical Chemistry B_ 2005, _109_ , 23955–23962, PMID: 16375384
  * Santiso and Trout 2011 Santiso, E. E.; Trout, B. L. A general set of order parameters for molecular crystals. _J. Chem. Phys._ 2011, _134_ , 064109
  * Shah et al. 2011 Shah, M.; Santiso, E. E.; Trout, B. L. Computer simulations of homogeneous nucleation of benzene from the melt. _The Journal of Physical Chemistry B_ 2011, _115_ , 10400–10412
  * Salvalaglio et al. 2012 Salvalaglio, M.; Vetter, T.; Giberti, F.; Mazzotti, M.; Parrinello, M. Uncovering molecular details of urea crystal growth in the presence of additives. _Journal of the American Chemical Society_ 2012, _134_ , 17221–17233
  * Salvalaglio et al. 2013 Salvalaglio, M.; Vetter, T.; Mazzotti, M.; Parrinello, M. Controlling and predicting crystal shapes: The case of urea. _Angew. Chem. Int._ 2013, _125_ , 13611–13614
  * Giberti et al. 2015 Giberti, F.; Salvalaglio, M.; Mazzotti, M.; Parrinello, M. Insight into the nucleation of urea crystals from the melt. _Chem. Eng. Sci._ 2015, _121_ , 51–59
  * Giberti et al. 2015 Giberti, F.; Salvalaglio, M.; Parrinello, M. Metadynamics studies of crystal nucleation. _IUCrJ_ 2015, _2_ , 256–266
  * Salvalaglio et al. 2015 Salvalaglio, M.; Perego, C.; Giberti, F.; Mazzotti, M.; Parrinello, M. Molecular-dynamics simulations of urea nucleation from aqueous solution. _Proc. Natl. Acad. Sci. U.S.A._ 2015, _112_ , E6–E14
  * Bjelobrk et al. 2019 Bjelobrk, Z.; Piaggi, P. M.; Weber, T.; Karmakar, T.; Mazzotti, M.; Parrinello, M. Naphthalene crystal shape prediction from molecular dynamics simulations. _CrystEngComm_ 2019, _21_ , 3280–3288
  * Bjelobrk et al. 2021 Bjelobrk, Z.; Mendels, D.; Karmakar, T.; Parrinello, M.; Mazzotti, M. Solubility prediction of organic molecules with molecular dynamics simulations. _Crystal Growth & Design_ 2021, _21_ , 5198–5205
  * Keys et al. 2011 Keys, A. S.; Iacovella, C. R.; Glotzer, S. C. Characterizing structure through shape matching and applications to self-assembly. _Annu. Rev. Condens. Matter Phys._ 2011, _2_ , 263–285
  * Shetty et al. 2002 Shetty, R.; Escobedo, F. A.; Choudhary, D.; Clancy, P. A novel algorithm for characterization of order in materials. _The Journal of chemical physics_ 2002, _117_ , 4000–4009
  * Duff and Peters 2011 Duff, N.; Peters, B. Polymorph specific RMSD local order parameters for molecular crystals and nuclei: $\alpha$-, $\beta$-, and $\gamma$-glycine. _The Journal of chemical physics_ 2011, _135_ , 134101
  * Nada 2020 Nada, H. Pathways for the formation of ice polymorphs from water predicted by a metadynamics method. _Scientific reports_ 2020, _10_ , 4708
  * Malkin et al. 2012 Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of ice crystallized from supercooled water. _Proceedings of the National Academy of Sciences_ 2012, _109_ , 1041–1045
  * Nettleton and Green 1958 Nettleton, R. E.; Green, M. S. Expression in Terms of Molecular Distribution Functions for the Entropy Density in an Infinite System. _The Journal of Chemical Physics_ 1958, _29_ , 1365–1370
  * Mendels et al. 2018 Mendels, D.; McCarty, J.; Piaggi, P. M.; Parrinello, M. Searching for Entropically Stabilized Phases: The Case of Silver Iodide. _Journal of Physical Chemistry C_ 2018, _122_ , 1786–1790
  * Piaggi and Parrinello 2018 Piaggi, P. M.; Parrinello, M. Predicting polymorphism in molecular crystals using orientational entropy. _Proceedings of the National Academy of Sciences of the United States of America_ 2018, _115_ , 10251–10256
  * Jones et al. 2001 Jones, E.; Oliphant, T.; Peterson, P. SciPy: Open Source Scientific Tools for Python. 2001,
  * Müllner 2013 Müllner, D. fastcluster: Fast Hierarchical, Agglomerative Clustering Routines for R and Python. _Journal of Statistical Software_ 2013, _53_ , 1–18
  * Cesa-Bianchi and Lugosi 2006 Cesa-Bianchi, N.; Lugosi, G. _Prediction, Learning, and Games_ ; Cambridge University Press, 2006
  * Amodeo et al. 2020 Amodeo, J.; Pietrucci, F.; Lam, J. Out-of-equilibrium polymorph selection in nanoparticle freezing. _Journal of Physical Chemistry Letters_ 2020, _11_ , 8060–8066
  * Gobbo et al. 2018 Gobbo, G.; Bellucci, M. A.; Tribello, G. A.; Ciccotti, G.; Trout, B. L. Nucleation of Molecular Crystals Driven by Relative Information Entropy. _Journal of Chemical Theory and Computation_ 2018, _14_ , 959–972
  * Lindley 1959 Lindley, D. V. Information Theory and Statistics. Solomon Kullback. New York: John Wiley and Sons, Inc. London: Chapman and Hall, Ltd. 1959. Pp. xvii, 395. $12.50. _Journal of the American Statistical Association_ 1959, _54_ , 825–827
  * Silverman 2018 Silverman, B. _Density Estimation for Statistics and Data Analysis_ ; Routledge, 2018
  * Song et al. 2020 Song, H.; Vogt-Maranto, L.; Wiscons, R.; Matzger, A. J.; Tuckerman, M. E. Generating Cocrystal Polymorphs with Information Entropy Driven by Molecular Dynamics-Based Enhanced Sampling. _Journal of Physical Chemistry Letters_ 2020, _11_ , 9751–9758
  * Invernizzi et al. 2017 Invernizzi, M.; Valsson, O.; Parrinello, M. Coarse graining from variationally enhanced sampling applied to the Ginzburg–Landau model. _Proceedings of the National Academy of Sciences_ 2017, _114_ , 3370–3374
  * Niu et al. 2018 Niu, H.; Piaggi, P. M.; Invernizzi, M.; Parrinello, M. Molecular dynamics simulations of liquid silica crystallization. _Proceedings of the National Academy of Sciences of the United States of America_ 2018, _115_ , 5348–5352
  * Lorch 1970 Lorch, E. Conventional and elastic neutron diffraction from vitreous silica. _Journal of Physics C: Solid State Physics_ 1970, _3_ , 1314–1322
  * Lin and Zhigilei 2006 Lin, Z.; Zhigilei, L. V. Time-resolved diffraction profiles and atomic dynamics in short-pulse laser-induced structural transformations: Molecular dynamics study. _Physical Review B_ 2006, _73_
  * Ramakrishnan and Yussouff 1979 Ramakrishnan, T. V.; Yussouff, M. First-principles order-parameter theory of freezing. _Physical Review B_ 1979, _19_ , 2775
  * Bonati and Parrinello 2018 Bonati, L.; Parrinello, M. Silicon Liquid Structure and Crystal Nucleation from Ab Initio Deep Metadynamics. _Physical Review Letters_ 2018, _121_ , 265701
  * Niu et al. 2019 Niu, H.; Yang, Y. I.; Parrinello, M. Temperature Dependence of Homogeneous Nucleation in Ice. _Physical Review Letters_ 2019, _122_ , 245501
  * Zhang et al. 2019 Zhang, Y.-Y.; Niu, H.; Piccini, G.; Mendels, D.; Parrinello, M. Improving collective variables: The case of crystallization. _J. Chem. Phys._ 2019, _150_ , 094509
  * Deng et al. 2021 Deng, Y.; Du, T.; Hui, L. Relationship of structure and mechanical property of silica with enhanced sampling and machine learning. _Journal of American Ceramic Society_ 2021, _104_ , 3910–3920
  * Ahlawat et al. 2021 Ahlawat, P.; Hinderhofer, A.; Alharbi, E. A.; Lu, H.; Ummadisingu, A.; Niu, H.; Invernizzi, M.; Zakeeruddin, S. M.; Dar, M. I.; Schreiber, F.; Hagfeldt, A.; Grätzel, M.; Rothlisberger, U.; Parrinello, M. A combined molecular dynamics and experimental study of two-step process enabling low-temperature formation of phase-pure $\alpha$FAPbI3. _Science Advances_ 2021, _7_ , eabe3326
  * Lodesani et al. 2021 Lodesani, F.; Tavanti, F.; Menziani, M. C.; Maeda, K.; Takato, Y.; Urata, S.; Pedone, A. Exploring the crystallization path of lithium disilicate through metadynamics simulations. _Phys. Rev. Materials_ 2021, _5_ , 075602
  * Lodesani et al. 2022 Lodesani, F.; Menziani, M. C.; Urata, S.; Pedone, A. Biasing crystallization in fused silica: An assessment of optimal metadynamics parameters. _Journal of Chemical Physics_ 2022, _156_ , 194501
  * Badin and Martoňák 2021 Badin, M.; Martoňák, R. Nucleating a Different Coordination in a Crystal under Pressure: A Study of the $B1\text{$-$}B2$ Transition in NaCl by Metadynamics. _Phys. Rev. Lett._ 2021, _127_ , 105701
  * Russo and Tanaka 2016 Russo, J.; Tanaka, H. Crystal nucleation as the ordering of multiple order parameters. _The Journal of Chemical Physics_ 2016, _145_ , 211801
  * Piccini et al. 2018 Piccini, G.; Mendels, D.; Parrinello, M. Metadynamics with discriminants: A tool for understanding chemistry. _Journal of chemical theory and computation_ 2018, _14_ , 5040–5044
  * Mendels et al. 2018 Mendels, D.; Piccini, G.; Brotzakis, Z. F.; Yang, Y. I.; Parrinello, M. Folding a small protein using harmonic linear discriminant analysis. _The Journal of chemical physics_ 2018, _149_ , 194113
  * Schwantes and Pande 2013 Schwantes, C. R.; Pande, V. S. Improvements in Markov state model construction reveal many non-native interactions in the folding of NTL9. _Journal of chemical theory and computation_ 2013, _9_ , 2000–2009
  * Pérez-Hernández et al. 2013 Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. Identification of slow molecular order parameters for Markov model construction. _The Journal of chemical physics_ 2013, _139_ , 07B604_1
  * M. Sultan and Pande 2017 M. Sultan, M.; Pande, V. S. tICA-metadynamics: accelerating metadynamics by using kinetically selected collective variables. _Journal of chemical theory and computation_ 2017, _13_ , 2440–2447
  * McCarty and Parrinello 2017 McCarty, J.; Parrinello, M. A variational conformational dynamics approach to the selection of collective variables in metadynamics. _The Journal of chemical physics_ 2017, _147_ , 204109
  * Zou et al. 2021 Zou, Z.; Tsai, S.-T.; Tiwary, P. Toward Automated Sampling of Polymorph Nucleation and Free Energies with the SGOOP and Metadynamics. _The Journal of Physical Chemistry B_ 2021, _125_ , 13049–13056
  * Tiwary and Berne 2016 Tiwary, P.; Berne, B. Spectral gap optimization of order parameters for sampling complex molecular systems. _Proceedings of the National Academy of Sciences_ 2016, _113_ , 2839–2844
  * Rogal et al. 2019 Rogal, J.; Schneider, E.; Tuckerman, M. E. Neural-Network-Based Path Collective Variables for Enhanced Sampling of Phase Transformations. _Phys. Rev. Lett._ 2019, _123_ , 245701
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. _Physical review letters_ 2007, _98_ , 146401
  * Behler 2011 Behler, J. Atom-centered symmetry functions for constructing high-dimensional neural network potentials. _The Journal of chemical physics_ 2011, _134_ , 074106
  * Bonati et al. 2020 Bonati, L.; Rizzi, V.; Parrinello, M. Data-driven collective variables for enhanced sampling. _J Phys Chem Lett_ 2020, _11_ , 2998–3004
  * Karmakar et al. 2021 Karmakar, T.; Invernizzi, M.; Rizzi, V.; Parrinello, M. Collective variables for the study of crystallisation. _Molecular Physics_ 2021, _119_ , e1893848
  * Bonati et al. 2021 Bonati, L.; Piccini, G.; Parrinello, M. Deep learning the slow modes for rare events sampling. _Proceedings of the National Academy of Sciences_ 2021, _118_ , e2113533118
  * Niu et al. 2019 Niu, H.; Yang, Y. I.; Parrinello, M. Temperature dependence of homogeneous nucleation in ice. _Phys. Rev. Lett._ 2019, _122_ , 245501
  * Sultan and Pande 2018 Sultan, M. M.; Pande, V. S. Automated design of collective variables using supervised machine learning. _The Journal of chemical physics_ 2018, _149_ , 094106
  * Chen et al. 2018 Chen, W.; Tan, A. R.; Ferguson, A. L. Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design. _The Journal of chemical physics_ 2018, _149_ , 072312