2024 AIChE Annual Meeting
(126e) Modeling Nucleation of Molecular Crystals from Solution
Nucleation and crystallization are ubiquitous in nature, and frequently occur in industrial settings. Some
familiar examples include formation of clouds, gout, and kidney stones. Industrial examples could be Portland
cement (the basic ingredient of concrete, stucco, and mortar), and a wide range of pharmaceuticals. It is
estimated over 90% of pharmaceuticals contain one or more crystalline substances.1
In drug manufacture, understanding nucleation and subsequent crystallization is essential to selectively
creating a uniform product. Nucleation is utilized at several stages in the manufacturing process, including
intermediate purification steps, or to control particle size. Several fundamental aspects of crystallization and
nucleation are poorly understood, despite their prevalence. Drug molecules can have a wide variety of func-
tional groups, and their interactions with each other, co-solutes, the solvent, anti-solvents, and impurities are
quite diverse. This large number of independent variables involved in a nucleation makes systematic exper-
imental study extremely challenging. For example, approximately 37% of single component compounds are
known to exhibit some degree of polymorphism, which is directly influenced by the initial nucleation event.2
Nucleation is a stochastic process, and can involve just a few ions.3, which compounds the aforementioned
experimental difficulties in the investigation of nucleation, and post-nucleation crystal growth. In spite of
these great challenges, recent progress has been made.4,5,6
Nucleation of solids from super-saturated solutions is a good candidate for molecular dynamics inves-
tigation which is the main focus of this work. Nucleation simulations are nevertheless a complex, difficult
problem. Typically, the timescale for nucleation is orders of magnitude greater than the timescales available
for simulation. In addition, it is difficult to define order parameters which can efficiently describe the tran-
sition between the two phases. Many popular biased techniques are limited to just a few collective variables
due to poor scaling for additional degrees of freedom added, such as umbrella sampling7, and metadynamics8.
Identifying a “path” collective variable which approximates the reaction coordinate via machine learning9,10
or Monte Carlo techniques is another option, but it requires significant effort and data for each new system.
In this work, we studied the nucleation of the model drug Sulfadiazine in a supersaturated solution with
acetonitrile. We looked at the free energy barrier to the formation of a critical nucleus at five different
degrees of supersaturation. In addition, we used classical nucleation theory to derive a heuristic point of
comparision for the energy barrier to nucleation, to validate our experimental results.
Methodology
Our biasing technique was the string method in collective variables, as defined by Maragliano et al11, and
a set of general molecular order parameters defined by Santiso and Trout.12 The string method only explores
local regions of phase space, not the entire collective variable space, resulting in excellent scaling for large
numbers of collective variables, as needed by the Santiso-Trout order parameters. In the string method, one
creates a series of replicas of the system along a path in collective varaible space, connecting two metastable
basins in phase space. Each replica is then allowed to relax, according to the negative gradient of the free
energy with respect to the collective variables at that point on the string. Periodically, a new string is created,
and new images are placed on it. This is to avoid having all the images ending up in the metastable basins.
The Santiso-Trout order parameters are able to completely describe different crystalline phases using an
"extended" pairwise distribution. This extended distribution incorporates additional directional probability
distributions beyond the typical gaussians, ultimately resulting in per-molecule order parameters which are
each sensitive to a different kind of symmetry.
Results
We studied the model drug sulfadiazine, at five different levels of supersaturation in acetonitrile. We
found that the free energy barrier to nucleation decreases monotonically with supersaturation, as does the
size of the critical nucleus. We drew from classical nucleation theory, and found a heuristic rule in terms
of saturation, and surface energy which was in good agreement with our simulations. Finally, we developed
an open source module implementing the aforementioned order parameters into PLUMED13 which, along
with our in-house suite of code, can be used for free energy calculations via the string method in collective
variables on any platform supporting PLUMED. As supersaturation decreased, we found instabilities, as the
cutoff radius of our order parameters was less than the average density of the solute molecules in solution.
We found that adding a local density order parameter greatly increased the stability of simulations, and
aided in convergence of the string.
Discussion
Understanding the trends in the energy barrier is a good first step to developing a more general procedure
for the study the nucleation of molecular crystals in solution. Just as moving forward from lennard-jones
particles to nucleation of molecules from the melt is a natural step, we have now moved to modeling nucleation
of a molecular crystal in solution. In the future this could help inform design choices for pharmaceutical
manufacture. Drugs are known to be extremely sensitive to the choice of solvent, and there is no widely
accepted rules for solvent choice. Understanding trends in the energy barrier to nucleation in different
solvents could be a first step in that process.
The plumed project is complete, and at time of writing there are no known bugs. This means our order
parameters are now available on a wide range of MD platforms, in a well documented manner. It is also well
parallelized, and far more expandable, and maintainable than previously, as all future modifications can be
done within the PLUMED ecosystem.
Conclusion
We explored the nucleation of the drug Sulfadizine at differing degrees of supersaturation, and found
plausible trends in that barrier, as validated by a classical model. We developed software to aid in future
works, and lower the barrier of entry for other groups to use our methods. We hope to continue our work
next by exploring the role of solvent in the nucleation, by applying our order parameters and the string
method to both the solute and solvent.
References
1. Shekunov, B.Y., York, Crystallization processes in pharmaceutical technology and drug delivery design.
Journal of Crystal Growth 211 (2000) 122}136 https://doi.org/10.1016/S0022-0248(99)00819-2
2. Saha BK, Nath NK, Thakuria R. Polymorphs with Remarkably Distinct Physical and/or Chemical
Properties. Chem Rec. 2023 Jan;23(1):e202200173. doi: 10.1002/tcr.202200173. Epub 2022 Sep 27.
PMID: 36166697.
3. Zimmermann, N.E.R. et al Nucleation of NaCl from Aqueous Solution: Critical Sizes, Ion-Attachment
Kinetics, and Rates J. Am. Chem. Soc.2015, 137, 41, 13352–13361 https://doi.org/10.1021/jacs.
5b08098
4. Liu, J. Direct Measurement of Critical Nucleus Size in Confined Volumes Langmuir 2007, 23, 13,
7286–7292 Publication Date:May 22, 2007 https://doi.org/10.1021/la063650a
5. Wang, Z., Wang, F., Peng, Y. et al. Direct observation of liquid nucleus growth in homogeneous
melting of colloidal crystals. Nat Commun 6, 6942 (2015). https://doi.org/10.1038/ncomms7942
6. Kalita, A., Mrozek-McCourt, M., Kaldawi, T.F. et al. Microstructure and crystal order during
freezing of supercooled water drops. Nature 620, 557–561 (2023). https://doi.org/10.1038/
s41586-023-06283-2
7. Kästner, J. Umbrella sampling WIREs Comput. Mol. Sci. 2011, 1, 932– 942 DOI: 10.1002/wcms.66
8. Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics WIREs Comput. Mol. Sci. 2011, 1, 826–
843 DOI: 10.1002/wcms.31
9. Ray, D. Deep Learning Collective Variables from Transition Path Ensemble. J. Chem. Phys. 158,
204102 (2023) https://doi.org/10.1063/5.0148872
10. Bonati, L. Data-Driven Collective Variables for Enhanced Sampling. Phys. Chem. Lett. 2020, 11, 8,
2998–300 https://doi.org/10.1021/acs.jpclett.0c00535
11. L. Maragliano et al., J. Chem. Phys. 125, 024106 (2006)
12. Santiso EE, Trout BL. A general set of order parameters for molecular crystals. J Chem Phys. 2011
Feb 14;134(6):064109. doi: 10.1063/1.3548889. PMID: 21322663.
13. G.A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, G. Bussi. PLUMED2: New feathers for an
old bird,