IntroductionNucleation and crystallization are ubiquitous in nature, and frequently occur in industrial settings. Somefamiliar examples include formation of clouds, gout, and kidney stones. Industrial examples could be Portlandcement (the basic ingredient of concrete, stucco, and mortar), and a wide range of pharmaceuticals. It isestimated over 90% of pharmaceuticals contain one or more crystalline substances.1In drug manufacture, understanding nucleation and subsequent crystallization is essential to selectivelycreating a uniform product. Nucleation is utilized at several stages in the manufacturing process, includingintermediate purification steps, or to control particle size. Several fundamental aspects of crystallization andnucleation 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 arequite 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 areknown to exhibit some degree of polymorphism, which is directly influenced by the initial nucleation event.2Nucleation is a stochastic process, and can involve just a few ions.3, which compounds the aforementionedexperimental difficulties in the investigation of nucleation, and post-nucleation crystal growth. In spite ofthese great challenges, recent progress has been made.4,5,6Nucleation 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, difficultproblem. Typically, the timescale for nucleation is orders of magnitude greater than the timescales availablefor 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 variablesdue 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,10or 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 withacetonitrile. We looked at the free energy barrier to the formation of a critical nucleus at five differentdegrees of supersaturation. In addition, we used classical nucleation theory to derive a heuristic point ofcomparision for the energy barrier to nucleation, to validate our experimental results.MethodologyOur biasing technique was the string method in collective variables, as defined by Maragliano et al11, anda set of general molecular order parameters defined by Santiso and Trout.12 The string method only exploreslocal regions of phase space, not the entire collective variable space, resulting in excellent scaling for largenumbers of collective variables, as needed by the Santiso-Trout order parameters. In the string method, onecreates a series of replicas of the system along a path in collective varaible space, connecting two metastablebasins in phase space. Each replica is then allowed to relax, according to the negative gradient of the freeenergy 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 probabilitydistributions beyond the typical gaussians, ultimately resulting in per-molecule order parameters which areeach sensitive to a different kind of symmetry.ResultsWe studied the model drug sulfadiazine, at five different levels of supersaturation in acetonitrile. Wefound that the free energy barrier to nucleation decreases monotonically with supersaturation, as does thesize of the critical nucleus. We drew from classical nucleation theory, and found a heuristic rule in termsof saturation, and surface energy which was in good agreement with our simulations. Finally, we developedan open source module implementing the aforementioned order parameters into PLUMED13 which, alongwith our in-house suite of code, can be used for free energy calculations via the string method in collectivevariables 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, andaided in convergence of the string.DiscussionUnderstanding the trends in the energy barrier is a good first step to developing a more general procedurefor the study the nucleation of molecular crystals in solution. Just as moving forward from lennard-jonesparticles to nucleation of molecules from the melt is a natural step, we have now moved to modeling nucleationof a molecular crystal in solution. In the future this could help inform design choices for pharmaceuticalmanufacture. Drugs are known to be extremely sensitive to the choice of solvent, and there is no widelyaccepted rules for solvent choice. Understanding trends in the energy barrier to nucleation in differentsolvents 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 orderparameters are now available on a wide range of MD platforms, in a well documented manner. It is also wellparallelized, and far more expandable, and maintainable than previously, as all future modifications can bedone within the PLUMED ecosystem.ConclusionWe explored the nucleation of the drug Sulfadizine at differing degrees of supersaturation, and foundplausible trends in that barrier, as validated by a classical model. We developed software to aid in futureworks, and lower the barrier of entry for other groups to use our methods. We hope to continue our worknext by exploring the role of solvent in the nucleation, by applying our order parameters and the stringmethod 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,