A complete description of thermodynamic stabilities of molecular crystals
AA complete description of thermodynamic stabilities of molecular crystals
Venkat Kapil ∗ Yusuf Hamied Department of Chemistry, University of Cambridge,Lensfield Road, Cambridge, CB2 1EW,UK andLaboratory of Computational Science and Modeling, Institut des Matériaux,École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
Edgar A Engel † TCM Group, Cavendish Laboratory, University of Cambridge,J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom (Dated: March 1, 2021)Accurate prediction of the stability of molecular crystals is a longstanding challenge, as oftenminuscule free energy differences between polymorphs are sensitively affected by the description ofelectronic structure, the statistical mechanics of the nuclei and the cell, and thermal expansion. Theimportance of these effects has been individually established, but rigorous free energy calculations,which simultaneously account for all terms, have been prevented by prohibitive computational costs.Here we reproduce the experimental stabilities of polymorphs of three prototypical compounds –benzene, glycine, and succinic acid – by computing rigorous ab initio
Gibbs free energies, at afraction of the cost of conventional harmonic approximations. This is achieved by a bottom-upapproach, which involves generating machine-learning potentials to calculate surrogate free energiesand subsequently calculating true ab initio free energies using inexpensive free energy perturbations.Accounting for all relevant physical effects is no longer a daunting task and provides the foundationfor reliable structure predictions for more complex systems of industrial importance.
Keywords: polymorphism, free energy, machine learning
Molecular crystals are ubiquitous in the phar-maceutical industry [1] and show great promise forapplications in organic photovoltaics [2], gas absorp-tion [3], and the food [4], pesticide [5] and fertilizerindustries [6]. Their tendency to exhibit polymorphismi.e. to exist in multiple crystal structures, on onehand provides a mechanism to tune properties bycontrolling crystal structure [7], and on the other handintroduces the challenge of synthesising and stabilizingcrystal structures with desired properties [8]. Whilethermodynamic stability at the temperature and pres-sure of interest is sufficient (although not necessary )to ensure long-term stability, simply understandingthermodynamic stability already poses a formidablechallenge. This is particularly true for pharmaceuticals,where free energy differences between drug polymorphsare often smaller than 1 kJ/mol [9], leading to therisk of the drug transforming into a less soluble andconsequently less effective form during the course ofmanufacturing, storage or shelf-life [10, 11]. Indeed,the problem of late appearing drug polymorphs iswidespread [12, 13].The pharmaceutical industry therefore spends con- ∗ Correspondence email address: [email protected] † Correspondence email address: [email protected] kinetics may protect thermodynamically metastable structuresfrom decaying almost indefinitely siderable resources on high-throughput crystallizationexperiments to screen for (meta-)stable, synthesizablepolymorphs [14], into which the target polymorphmay decay. However, crystallization experiments donot probe thermodynamic stability, and conclusivestudies of the impact of temperature changes aftercrystallisation on the stability of polymorphs (i.e. theirmonotropic or enantiotropic nature [15]) are often pre-vented by limited sample quantities. Hence the appealof theoretical crystal structure prediction (CSP) [16]based on the thermodynamic stability of polymorphs,which promises to guide industrial drug design and tocomplement crystallization experiments [17].Despite the demonstrable value of CSP for manyclasses of materials [21–32], and the continuing progressevidenced by a series of blind tests [33–38], the successof CSP for molecular crystals has been limited bythe inability to routinely predict the relative stabil-ity of competing candidate structures [38]. This islargely because the methods used for stability rankingstypically ignore or approximate the subtle interplayof several effects, such as intricate inter-molecularinteractions [39], the (quantum) statistical mechanicsof the nuclei [40] and the unit cell [41], and thermalexpansion [42], thereby incurring errors larger thanthe free energy differences of interest. Indeed, theimportance of these effects has been demonstrated,albeit in isolation. For instance, harmonic [43–46] andanharmonic [40, 41] vibrational free energies have been a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b Figure 1. Structures of forms I and II of benzene containing16 molecules, forms α , β , and γ of glycine containing 24molecules, and forms α and β of succinic acid containing 24molecules. Hydrogen, carbon, nitrogen, and oxygen atomsare shown in white, gray, blue, and red, respectively. studied on the basis of force-fields or semi-local DFT,thereby ignoring their interplay with other importanteffects like Pauli repulsion, many-body dispersion,and fluctuations of the simulation cell. While all ofthese effects can be simultaneously and accurately beaccounted for using ab initio path integral (PI) freeenergy calculations [47], the cost of such simulationshas generally necessitated the use of less accurate butcomputationally cheaper approximate methods [48].This work addresses the problem of accurately rank-ing a preexisting set of polymorphs of complex molec-ular crystals by thermodynamic stability. Following inthe footsteps of Ref. [49], we calculate accurate ab initio Gibbs free energies of diverse, complex molecular crys-tals while leveraging accurate machine-learning poten-tials [50] (MLPs) to mitigate computational costs. Wedemonstrate that computing MLP free energies whileexactly accounting for the statistical mechanics of thenuclei and the cell, and subsequently promoting them to
Figure 2. Path integral (PI) Gibbs free energy differencesbetween forms I and II of benzene (B-I and B-II), α , β , and γ -glycine (G- α , G- β and G- γ ), and α and β -succinic acid(S- α and S- β ) calculated using PBE0-MBD based machine-learning potentials (MLPs) (blue) and corrected to the abinitio PBE0-MBD DFT level using free energy perturbation(red). Experimental data [18–20] are shown in green. ab initio hybrid-functional DFT free energies using freeenergy perturbation (FEP), leads to consistently cor-rect stability predictions for polymorphs of three proto-typical molecular crystals – benzene, glycine, and suc-cinic acid – and thus constitutes an avenue to predic-tive CSP for complex molecular crystals of industrialimportance. While we focus on ambient pressure andtemperatures of 100 K and 300 K, respectively, we alsocompute gradients of Gibbs free energies as indicatorsfor the monotropic or enantiotropic nature of the poly-morphs, and note that our free energy calculations nat-urally generalize to the calculation of full p − T phasediagrams. RESULTS
In order to predict rigorous relative stabilities, wecombine PI thermodynamic integration [40, 48] in theconstant pressure ensemble (thereby accounting for an-harmonic quantum nuclear motion and the fluctuationsand thermal expansion of the cell) with DFT calcula-tions with the hybrid PBE0 functional and the many-body dispersion correction of Tkatchenko et al. (inthe following referred to as PBE0-MBD). PBE0-MBD
Figure 3. Schematic representation of the workflow for computing ab initio , quantum anharmonic Gibbs free energies forcandidate crystal structures. The upper half of the figure shows the main steps: (1) generating ab initio reference data onwhich to (2) train a combined MLP, which can then be used to (3) compute MLP Gibbs free energies, which one can finally(4) promote to ab initio
Gibbs free energies. The lower half (shaded in blue) details the key aspects of how each of thesesteps is performed in practice. provides an accurate description of inter-molecular in-teractions, as benchmarked using experimental andCCSD(T) lattice energies for various molecular crys-tals, including form I of benzene and α -glycine [51, 52].Since direct calculation of Gibbs free energies from ab initio PI thermodynamic integration is preventedby the cost of the required energy and force evalua-tions [48], ab initio
Gibbs free energies are calculatedin a four-step process, as depicted schematically inFig. 3. The first and second step involve the generationof PBE0-MBD reference data and the construction ofMLPs thereupon. Both are detailed in the methodssection, use readily available and well-documentedsoftware, and can easily be repeated for almost anycompound of interest at the cost of a few thousandPBE0-MBD unit-cell reference calculations. Crucially,accurate MLPs, which predict energies and forces ofextended structures on the basis of local contribu-tions, do not necessarily require reference data forthe extended structures, provided all relevant localenvironments are included [53]. In cases which do notbenefit from the small unit cell sizes considered in thiswork, reference data for smaller “disordered” structures should thus provide an alternative to explicit ab initio calculations for large crystal structures. In a thirdstep, Gibbs free energies for the much larger simulationsupercells are calculated using PI thermodynamicintegrations based on the surrogate MLPs, which aremany orders of magnitude cheaper to evaluate than the ab initio reference method. However, residual errorswith respect to the ab initio free energies arise fromthe imperfect reproduction of the ab initio potentialenergy surface by the MLPs [49, 54, 55]. The typicalerrors in MLP predictions of configurational energiesare shown in Table I. These may arise from the lack oflong-range interactions in the parametrization of theMLP [56], the incomplete description of local atomicenvironments [57], or simply due to insufficient trainingdata or the stochastic nature of the training procedure.In a fourth and final step we eliminate the associatederrors to obtain true ab initio
Gibbs free energiesby computing the difference between the MLP andPBE0-MBD free energies using FEP [49].We demonstrate the effectiveness of the scheme ona diverse set of prototypical systems: benzene is thearchetypal rigid, van-der-Waal’s bonded molecular
System reference data energy RMSE [kJ/mol]Benzene 1,000 1.2Glycine 4,000 1.6Succinic acid 2,000 2.3Table I. Number of single-point PBE0-MBD calculations un-derlying each MLP, and their respective root-mean-squareerrors (RMSE) in predicting energies on a separate test setof configurations from PI simulations of the experimentalunit cells. crystal, while glycine prototypes flexible zwitter-ionicsystems, and succinic acid represents general hydrogen-bonded systems. We focus on free energy differences atambient pressure and 100 K for benzene and succinicacid and 300 K for glycine, but full p - T phase-diagramscan be calculated analogously. For each compoundwe consider the stable polymorph and its closestcompetitors: forms I and II of benzene [58], α , β , and γ -glycine [59], and α and β -succinic acid [60]. Thenearly orthorhombic simulation supercells shown inFig. 1, which contain equivalent numbers of moleculesfor all polymorphs of the same compound, ensureuniform sampling of the vibrational Brillouin zone(BZ), and that the centre-of-mass free energy termlargely cancels out in free energy differences betweenpolymorphs.As shown in Fig. 2, the final ab initio Gibbs freeenergies (shown in black) reproduce the greater sta-bility of form I over form II of benzene [61], of α over β -succinic acid [62], the metastability of β -glycine [20],and the near degeneracy of α and γ -glycine [19].Moreover, our Gibbs free energy differences are in goodagreement with calorimetry data for glycine to withinstatistical and experimental uncertainties. Meanwhile,the MLP-based stability predictions (shown in red) areonly limited by the accuracy with which the MLPsreproduce the ab initio potential energy surface (seeTable I) and, consequently, correctly reproduce thegreater stability of form I of benzene compared toform II. At the same time, the incorrect MLP-basedstability predictions for succinic acid and glycinehighlight the critical importance of the final FEPstep. Glycine and succinic acid are zwitter-ionic andhave delocalized ions, respectively, necessitating theexplicit treatment of long-ranged interactions. This isnot provided by the “local” MLP framework employedhere. Note that promoting MLP free energies tothe ab initio level by FEP only incurs the cost ofa few tens of ab initio energy and force evaluationsfor configurations sampled by the MLPs, and thusconstitutes a reliable and computationally efficientmeans of predicting the relative stability of polymorphs. Figure 4. MLP (free) energy differences between forms I andII of benzene (B-I and B-II), α , β , and γ -glycine (G- α , G- β and G- γ ), and α and β -succinic acid (S- α and S- β ) at differ-ent tiers of accuracy: PBE-TS level static-lattice energy dif-ference (Umin) for the optimized cells (green), bare PBE-TSbased QHA (gray) and single-point PBE0-MBD correctedQHA (pink), full PBE0-MBD based QHA (brown), and theexact PI free energy difference (blue). The shaded regionindicates free energy differences within 1 kJ/mol of the re-spective exact PI result as a guide to the eye. In order to highlight the advantages of choosing theapproach proposed here over established and stream-lined approximate methods for ranking stabilities inCSP, we assess the limitations of the most widely usedapproximate methods, prefaced by acknowledging theirsuccesses for a wide range of applications [39, 43, 45].To this end we restrict ourselves to MLP free energiesand compare the respective approximate (free) energydifferences between polymorphs to the correspondingexact MLP Gibbs free energies. The most widespreadapproach is to use lattice energies of variable-cellgeometry-optimized structures as a proxy for free en-ergies [38]. These are typically calculated using semi-local DFT, for instance with the PBE functional anda TS dispersion-correction (PBE-TS) [62–64]. Fig. 4shows that this incurs errors in free energy differencesup to above 2 kJ/mol. A second, more accurate ap-proach is to estimate finite-temperature free energieswithin the quasi-harmonic approximation (QHA) [65].For our systems this reduces the typical errors to justabove 1 kJ/mol (see Fig. 4). The arguably most reli-able strategy employed to date [43] is to correct QHAfree energies obtained using semi-local DFT based on asingle, hybrid-functional DFT calculation [39, 43]. Forour systems this correction only leads to a small im-provement, with the typical error not dipping signifi-cantly below 1 kJ/mol. Unfortunately, none of theseapproaches reproduce the correct stability order for allsystems. Given that errors of 1 kJ/mol are often consid-ered to be within “chemical accuracy”, it is worth em-phasising that the compounds considered here are nothand-picked, “pathological” examples, but expected tobe representative of many (bio-)molecular compounds.Notably, performing the full QHA at the hybrid-functional level (to account for the interplay of Paulirepulsion and many-body electronic dispersion withnuclear quantum effects and thermal expansion withina harmonic approximation) halves the average errorin free energy differences, and seemingly makes theQHA competitive with the rigorous PI approach. Itis therefore worthwhile to put the respective costs ofthe calculations into perspective. For glycine, as themost costly example, the 4,000 unit-cell PBE0-MBDcalculations constituting the reference data for theMLP add to the MLP-based PI thermodynamic inte-gration, and the 50 supercell PBE0-MBD calculationsrequired for the FEP to result in a total cost of around148,000 core-hours per polymorph. For comparison, anaive harmonic approximation using finite differencesat the same level of theory and for the same simulationsupercell would require 1,440 supercell PBE0-MBDcalculations costing around 1.7 million core-hours perpolymorph (which can be reduced by around a factorof four by using non-diagonal supercells to probeindividual k -points [66]). Despite a focus on universalapplicability over efficiency, the cost the above rigorousGibbs free energies is thus small compared to theestimated cost of calculating PBE0-MBD free energieseven within the harmonic approximation.Although polymorph stability at ambient conditionsis of particular significance, the phase behaviour at gen-eral thermodynamic conditions is of great interest forcontrolling polymorph crystallization. Temperature ar-guably constitutes the most easily tuneable thermody-namic constraint and is often key to crystallizing differ-ent polymorphs from the melt [15]. While we reserveevaluations of the full temperature phase-diagrams forfuture work, gradients of Gibbs free energies, includ-ing the molar volume, entropy, and constant-pressureheat capacity, come as a compliment of the above cal-culations and provide indication regarding prospectivechanges in relative stability and thus the monotropic orenantiotropic nature of polymorphs. Indeed these quan-tities are directly related to several “thermodynamicrules of thumb” for determining stability trends from ex-perimental densities, melting temperatures, and heats of transformation and fusion [15]. For instance, approx-imating Gibbs free energy differences to depend linearlyon pressure, ∆ G ( p, T ) ≈ ∆ G ( p , T ) + ( p − p )∆ V , wecan estimate form II of benzene to become thermody-namically over the ambient pressure form I at 1.4 GPa(at 100 K), which is in good agreement with the experi-mentally determined transition pressure of 1.5 GPa [61].Similarly, we determine the entropies of γ -glycine and α -succinic acid to be smaller than those of α -glycine and β -succinic acid, respectively, making the latter the pre-ferred high-temperature polymorphs in agreement withthe experimentally observed transitions to α -glycineand β -succinic acid at 443 K [67] and 410 K [62], re-spectively. DISCUSSION
Marrying state-of-the-art electronic structure, freeenergy, and machine-learning methods in a widely appli-cable framework renders rigorous, predictive free energycalculations not just computationally viable, but effi-cient. On top of highly accurate Gibbs free energies atthe chosen temperatures and pressures, our simulationsprovide (for free) entropies and molar volumes, whichconstitute valuable indicators of prospective phase-transitions, in direct analogy with the “thermodynamicrules of thumbs” used by experimentalists to assess themonotropic or enantiotropic nature of polymorphs. Fur-ther, a rigorous assessment of thermally-induced phase-transitions can be obtained at the cost of an additionaltemperature thermodynamic integration for each poly-morph. Indeed, our framework naturally generalizes tothe calculation of full p − T phase diagrams.Notably, the resultant predictions of thermodynamicstability are very robust with respect to the natureof the candidate polymorphs, since dynamic disorder,thermal expansion, conformational relaxation of themolecular units within the crystal structure, and poten-tial (dynamic) instabilities of candidate polymorphs areautomatically accounted for. For instance, our calcula-tions account for the suppressed rotation of the aminegroups in glycine [68] as well as the dynamic disorder ofthe protons in succinic acid. Furthermore, the dynam-ical nature of our simulations collapses redundant setsof minimum potential energy structures, as are com-mon in CSP, onto the typically much smaller number ofassociated free energy minima [69] (prior to full freeenergy computation), thereby eliminating a commonheadache in traditional CSP studies. Finally, calcu-lating Gibbs free energies at experimental conditionsallows direct comparison with experimental data, elim-inating the need to extrapolate experimental data to 0 Kand subtracting zero-point effects in an ad-hoc fashion,which is prone to errors.Despite the demonstrated power of the describedscheme for free energy calculations, some limitationsto its universality must be acknowledged, starting withthe empiricism involved in selecting the exchange-correlation functional and dispersion correction usedin the DFT calculations. Indeed, truly resolving thegreater stability of γ -glycine compared to α -glycine mayrequire a more accurate description of electronic struc-ture. The framework extends naturally to predictionsof Gibbs free energies based on quantum-chemical elec-tronic structure methods such as MP2 [70], RPA [71],coupled cluster [72] or quantum Monte Carlo [73], someof which are systematically improvable, and can therebybe rendered truly ab initio . While these come at an in-creased computational cost per calculation, ∆ -learningschemes [74] combined with active learning strate-gies [75], promise to minimize the number of quantum-chemical calculations required to train accurate MLPsand thus to keep the overall costs in check. Mean-while, uncertainty estimation for MLP predictions [76]may facilitate determining if a FEP correction is re-quired. Recent investigations of the incompleteness ofrepresentations of atomic structures [57], the use ofhigher body-order corrections [77], and means of includ-ing long-ranged interactions [56] pave the way for MLPscapable of reproducing quantum-chemical energies withsub-kJ/mol accuracy, promising to eliminate the needfor FEP to the reference level of theory.Indeed, the observation that MLP-based free ener-gies already reproduce the reference free energy differ-ences to within around 1 kJ/mol suggests that MLPsare already suitable for screening CSP data at the levelof configurational energies, as recently demonstrated inRef. [78], as well as more accurate free energies, priorto predictive ab initio free energy calculations for themost promising candidates. We thus envisage the use ofMLP-based free energy calculations to become commonpractice in CSP applications.The demonstrated framework promises to recover ac-curate stability rankings of more flexible and more chal-lenging systems such as salts, hydrates and co-crystalsof pharmaceutical importance, and functional materialssuch as metal organic frameworks, and porous organiccrystals. More insights on how a particular molecularcrystal can be synthesized, or quantification of decayrates into more stable structures can be obtained bystudying kinetic stability. The demonstrated accuracyof MLP-based free energy calculations sets the stage forfuture studies of kinetic effects as well as full p - T phase-diagrams using MLPs in a reliable and computationallyefficient manner. ACKNOWLEDGEMENTS
VK acknowledges funding from the Swiss NationalScience Foundation (SNSF), Project P2ELP2_191678, and support from the NCCR MARVEL, funded by theSNSF. EAE acknowledges funding from Trinity College,Cambridge. The authors thank Benjamin Shi, DaanFrenkel, Angelos Michaelides, Sally Price and MicheleCeriotti for valuable suggestions on the manuscript.EAE and VK acknowledges allocation of CPU hoursby CSCS under Project IDs s960 and s1000.
METHODS
Machine-learning potentials:
We have constructed Behler-Parinello type neural network potentials [79] for benzene,glycine, and succinic acid using the n2p2 code [80]. In thisframework, structures are encoded in terms of local atom-centered symmetry functions (SF) [79]. Initial sets of SFswere generated following the recipe of Ref. [81]. Based on thesame reference structure-property data subsequently used fortraining, the 128 (benzene and succinic acid) and 256 (glycine)most informative SFs were extracted via PCovCUR selection [82].Our data is based on Langevin-thermostatted PI NVT simula-tions at 300 K, performed using the i-Pi force engine [83] coupledto DFTB+ [84] calculations with the 3ob parametrization [85].For each polymorph multiple cells were simulated, rescalingthe experimental cell lengths and angles by up to 10% and5%, respectively. All trajectories for a given compound wereconcatenated and farthest-point sampled [86–88] to extractthe most distinct configurations for feature selection and MLPtraining. Subsequently, ab initio reference energies and forceswere evaluated for said configurations.To minimise the computational cost of the reference calculationsthe MLPs are composed of a baseline potential trained to re-produce energies and forces from more affordable PBE-DFT [89]calculations with a TS dispersion correction [90] (PBE-TS), anda ∆ -learning [74] correction trained (on ten times fewer trainingdata) to reproduce the difference between the baseline and moreexpensive calculations with the hybrid PBE0 functional [91, 92]and the MBD dispersion correction [93, 94] (PBE0-MBD). For aseparate test set, the MLPs reproduce PBE0-MBD energies withroot-mean-square errors of 1.2 kJ/mol for benzene, 1.6 kJ/molfor glycine, and 2.3 kJ/mol for succinic acid, respectively. A b initio DFT calculations: BE0+MBD calculations wereperformed using FHI-aims [95–97] with the standard FHI-aims“intermediate” basis sets and a Monkhorst-Pack k-point grid [98]with a maximum spacing of . × π Å − . The PBE-TSbaseline calculations for a ∆ -learning approach were performedusing Quantum Espresso v6.3, the same k-point grid, a wave-function cut-off energy of 100 Rydberg, and the optimised,norm-conserving Vanderbilt pseudopotentials from Ref. [99]. Free energy methods:
For each polymorph the average cellwas determined using MLP based path integral (PI) NST simu-lations [100] at the desired inverse temperature β , accounting foranharmonic quantum nuclear motion and anisotropic cell fluc-tuations. The difference between the Gibbs and Helmholtz freeenergies is computed from a MLP based PI NPT simulation basedon its average cell is G MLP ( P ext , β ) − A MLP ( V, β ) = P ext V + β − ln ρ ( V | P ext , β ) , where ρ ( V | P ext , β ) is the probability of observing the cell vol-ume V at external pressure P ext and inverse temperature β .A standard Kirkwood construction [101] that transforms the Hamiltonian from a harmonic to an anharmonic one provides thedifference between the anharmonic and the harmonic quantumHelmholtz free energies: A MLP ( V, β ) − A harMLP ( V, β ) = (cid:90) dλ (cid:68) ˆ H MLP − ˆ H harMLP (cid:69) V,β, ˆ H λ , where ˆ H λ is the Hamiltonian of the MLP alchemical system withthe potential U λ ≡ λU MLP + (1 − λ ) U harMLP , and (cid:104)·(cid:105) is the ensem-ble average computed from a PI NVT simulation. The referenceabsolute harmonic Helmholtz free energy is obtained from a har-monic approximation using A harMLP ( V, β ) = U MLP ( V ) + (cid:88) i (cid:20) (cid:126) ω i + β − ln (cid:16) − e − β (cid:126) ω i (cid:17)(cid:21) , where ω i is the frequency of the i -th phonon mode. In a finalstep, the ab initio Gibbs free energy is obtained from its MLPcounterpart by free energy perturbation using G ( P ext , β ) − G MLP ( P ext , β ) = − β − ln (cid:68) e − β ( U − U MLP ) (cid:69) P ext ,β, ˆ H MLP . [1] S. Datta and D. J. W. Grant, Nature Reviews DrugDiscovery , 42 (2004).[2] S. R. Forrest, Nature , 911 (2004).[3] T. Tozawa, J. T. A. Jones, S. I. Swamy, S. Jiang,D. J. Adams, S. Shakespeare, R. Clowes, D. Brad-shaw, T. Hasell, S. Y. Chong, C. Tang, S. Thomp-son, J. Parker, A. Trewin, J. Bacsa, A. M. Z. Slawin,A. Steiner, and A. I. Cooper, Nature Materials , 973(2009).[4] H. Hondoh and S. Ueno, Progress in Crystal Growthand Characterization of Materials Special Issue: Re-cent Progress on Fundamentals and Applications ofCrystal Growth; Proceedings of the 16th InternationalSummer School on Crystal Growth (ISSCG-16), ,398 (2016).[5] J. Yang, C. T. Hu, X. Zhu, Q. Zhu, M. D. Ward, andB. Kahr, Angewandte Chemie , 10299 (2017).[6] K. Honer, E. Kalfaoglu, C. Pico, J. McCann, andJ. Baltrusaitis, ACS Sustainable Chemistry & Engi-neering , 8546 (2017).[7] D. Gentili, M. Gazzano, M. Melucci, D. Jones, andM. Cavallini, Chemical Society Reviews , 2502(2019).[8] D. Chistyakov and G. Sergeev, Pharmaceutics (2020), 10.3390/pharmaceutics12010034.[9] A. J. Cruz-Cabeza, S. M. Reutzel-Edens, and J. Bern-stein, Chemical Society Reviews , 8619 (2015).[10] S. R. Chemburkar, J. Bauer, K. Deming, H. Spiwek,K. Patel, J. Morris, R. Henry, S. Spanton, W. Dziki,W. Porter, J. Quick, P. Bauer, J. Donaubauer, B. A.Narayanan, M. Soldani, D. Riley, and K. McFar-land, Organic Process Research & Development , 413(2000).[11] K. R. Chaudhuri, Expert Opinion on Drug Delivery ,1169 (2008).[12] I. B. Rietveld and R. Céolin, Journal of Pharmaceuti-cal Sciences , 4117 (2015).[13] D.-K. Bučar, R. W. Lancaster, and J. Bernstein,Angewandte Chemie International Edition , 6972(2015).[14] S. L. Morissette, O. Almarsson, M. L. Peterson, J. F.Remenar, M. J. Read, A. V. Lemmo, S. Ellis, M. J.Cima, and C. R. Gardner, Advanced Drug DeliveryReviews , 275 (2004). [15] E. H. Lee, Asian Journal of Pharmaceutical Sciences , 163 (2014).[16] S. L. Price, Chemical Society Reviews , 2098 (2014).[17] J. Nyman and S. M. Reutzel-Edens, Faraday Discus-sions , 459 (2018).[18] G. L. Perlovich, L. K. Hansen, and A. Bauer-Brandl,Journal of Thermal Analysis and Calorimetry , 699(2001).[19] V. A. Drebushchak, Y. A. Kovalevskaya, I. E. Paukov,and E. V. Boldyreva, Journal of Thermal Analysis andCalorimetry , 109 (2003).[20] V. A. Drebushchak, E. V. Boldyreva, Y. A. Ko-valevskaya, I. E. Paukov, and T. N. Drebushchak,Journal of Thermal Analysis and Calorimetry , 65(2005).[21] C. J. Pickard and . J. Needs, Nature Physics , 473(2007).[22] B. Monserrat, R. J. Needs, E. Gregoryanz, and C. J.Pickard, Phys. Rev. B , 134101 (2016).[23] Y. M. Ma, M. Eremets, A. R. Oganov, Y. Xie, I. Tro-jan, S. Medvedev, A. O. Lyakhov, M. Valle, andV. Prakapenka, Nature , 182 (2009).[24] Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, J. Chem.Phys. , 174712 (2014).[25] I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J.Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri,Phys. Rev. Lett. , 157004 (2015).[26] A. P. Drozdov, P. P. Kong, V. S. Minkov, S. P. Besedin,M. A. Kuzovnikov, S. Mozaffari, L. Balicas, F. F. Bal-akirev, D. E. Graf, V. B. Prakapenka, E. Greenberg,D. A. Knyazev, M. Tkacz, and M. I. Eremets, Nature , 528 (2019).[27] H. L. Zhuang and R. G. Hennig, Journal of the Min-erals, Metals & Materials Society , 366 (2014).[28] N. Mounet, M. Gibertini, P. Schwaller, D. Campi,A. Merkys, A. Marrazzo, T. Sohier, I. Eligio Castelli,A. Cepellotti, G. Pizzi, and N. Marzari, Nature Nan-otechnology , 246 (2018).[29] A. Marrazzo, M. Gibertini, D. Campi, N. Mounet, andN. Marzari, Phys. Rev. Lett. , 117701 (2018).[30] P. E. A. Turchi, P. Söderlind, and A. I. Landa, Journalof the Minerals, Metals & Materials Society , 375(2014).[31] B. D. Conduit, N. G. Jones, H. J. Stone, and G. J.Conduit, Materials & Design , 358 (2017). [32] P. Sarker, T. Harrington, C. Toher, C. Oses,M. Samiee, J.-P. Maria, D. W. Brenner, K. S. Vec-chio, and S. Curtarolo, Nature Communications ,4980 (2018).[33] J. P. M. Lommerse, W. D. S. Motherwell, H. L. Am-mon, J. D. Dunitz, A. Gavezzotti, D. W. M. Hof-mann, F. J. J. Leusen, W. T. M. Mooij, S. L. Price,B. Schweizer, M. U. Schmidt, B. P. v. Eijck, P. Verwer,and D. E. Williams, Acta Crystallographica Section B:Structural Science , 697 (2000).[34] W. D. S. Motherwell, H. L. Ammon, J. D. Dunitz,A. Dzyabchenko, P. Erk, A. Gavezzotti, D. W. M. Hof-mann, F. J. J. Leusen, J. P. M. Lommerse, W. T. M.Mooij, S. L. Price, H. Scheraga, B. Schweizer, M. U.Schmidt, B. P. v. Eijck, P. Verwer, and D. E. Williams,Acta Crystallographica Section B: Structural Science , 647 (2002).[35] G. M. Day, W. D. S. Motherwell, H. L. Ammon,S. X. M. Boerrigter, R. G. Della Valle, E. Venuti,A. Dzyabchenko, J. D. Dunitz, B. Schweizer, B. P. vanEijck, P. Erk, J. C. Facelli, V. E. Bazterra, M. B. Fer-raro, D. W. M. Hofmann, F. J. J. Leusen, C. Liang,C. C. Pantelides, P. G. Karamertzanis, S. L. Price,T. C. Lewis, H. Nowell, A. Torrisi, H. A. Scheraga,Y. A. Arnautova, M. U. Schmidt, and P. Verwer,Acta Crystallographica Section B: Structural Science , 511 (2005).[36] G. M. Day, T. G. Cooper, A. J. Cruz-Cabeza, K. E.Hejczyk, H. L. Ammon, S. X. M. Boerrigter, J. S. Tan,R. G. Della Valle, E. Venuti, J. Jose, S. R. Gadre,G. R. Desiraju, T. S. Thakur, B. P. van Eijck, J. C.Facelli, V. E. Bazterra, M. B. Ferraro, D. W. M. Hof-mann, M. A. Neumann, F. J. J. Leusen, J. Kendrick,S. L. Price, A. J. Misquitta, P. G. Karamertzanis,G. W. A. Welch, H. A. Scheraga, Y. A. Arnautova,M. U. Schmidt, J. van de Streek, A. K. Wolf, andB. Schweizer, Acta Crystallographica Section B: Struc-tural Science , 107 (2009).[37] D. A. Bardwell, C. S. Adjiman, Y. A. Arnautova,E. Bartashevich, S. X. M. Boerrigter, D. E. Braun,A. J. Cruz-Cabeza, G. M. Day, R. G. Della Valle,G. R. Desiraju, B. P. van Eijck, J. C. Facelli, M. B.Ferraro, D. Grillo, M. Habgood, D. W. M. Hofmann,F. Hofmann, K. V. J. Jose, P. G. Karamertzanis, A. V.Kazantsev, J. Kendrick, L. N. Kuleshova, F. J. J.Leusen, A. V. Maleev, A. J. Misquitta, S. Mohamed,R. J. Needs, M. A. Neumann, D. Nikylov, A. M.Orendt, R. Pal, C. C. Pantelides, C. J. Pickard, L. S.Price, S. L. Price, H. A. Scheraga, J. van de Streek,T. S. Thakur, S. Tiwari, E. Venuti, and I. K. Zhitkov,Acta Crystallographica Section B: Structural Science , 535 (2011).[38] A. M. Reilly, R. I. Cooper, C. S. Adjiman, S. Bhat-tacharya, A. D. Boese, J. G. Brandenburg, P. J. By-grave, R. Bylsma, J. E. Campbell, R. Car, D. H.Case, R. Chadha, J. C. Cole, K. Cosburn, H. M.Cuppen, F. Curtis, G. M. Day, R. A. DiStasio Jr,A. Dzyabchenko, B. P. van Eijck, D. M. Elking, J. A.van den Ende, J. C. Facelli, M. B. Ferraro, L. Fusti-Molnar, C.-A. Gatsiou, T. S. Gee, R. de Gelder,L. M. Ghiringhelli, H. Goto, S. Grimme, R. Guo, D. W. M. Hofmann, J. Hoja, R. K. Hylton, L. Iuz-zolino, W. Jankiewicz, D. T. de Jong, J. Kendrick,N. J. J. de Klerk, H.-Y. Ko, L. N. Kuleshova, X. Li,S. Lohani, F. J. J. Leusen, A. M. Lund, J. Lv,Y. Ma, N. Marom, A. E. Masunov, P. McCabe, D. P.McMahon, H. Meekes, M. P. Metz, A. J. Misquitta,S. Mohamed, B. Monserrat, R. J. Needs, M. A. Neu-mann, J. Nyman, S. Obata, H. Oberhofer, A. R.Oganov, A. M. Orendt, G. I. Pagola, C. C. Pan-telides, C. J. Pickard, R. Podeszwa, L. S. Price, S. L.Price, A. Pulido, M. G. Read, K. Reuter, E. Schnei-der, C. Schober, G. P. Shields, P. Singh, I. J. Sug-den, K. Szalewicz, C. R. Taylor, A. Tkatchenko, M. E.Tuckerman, F. Vacarro, M. Vasileiadis, A. Vazquez-Mayagoitia, L. Vogt, Y. Wang, R. E. Watson, G. A.de Wijs, J. Yang, Q. Zhu, and C. R. Groom, ActaCrystallographica Section B: Structural Science, Crys-tal Engineering and Materials , 439 (2016).[39] N. Marom, R. A. DiStasio, V. Atalla, S. Levchenko,A. M. Reilly, J. R. Chelikowsky, L. Leiserowitz,and A. Tkatchenko, Angewandte Chemie InternationalEdition , 6629 (2013).[40] M. Rossi, P. Gasparotto, and M. Ceriotti, PhysicalReview Letters , 115702 (2016).[41] E. Schneider, L. Vogt, and M. E. Tuckerman, ActaCrystallographica Section B: Structural Science, Crys-tal Engineering and Materials , 542 (2016).[42] H.-Y. Ko, R. A. DiStasio, B. Santra, and R. Car,Physical Review Materials , 055603 (2018).[43] J. Hoja, H.-Y. Ko, M. A. Neumann, R. Car, R. A.DiStasio, and A. Tkatchenko, Science Advances ,eaau3338 (2019).[44] B. P. van Eijck, W. T. M. Mooij, and J. Kroon, TheJournal of Physical Chemistry B , 10573 (2001).[45] A. M. Reilly and A. Tkatchenko, Physical Review Let-ters , 055701 (2014).[46] J. Nyman and G. M. Day, CrystEngComm , 5154(2015).[47] D. Marx and M. Parrinello, The Journal of ChemicalPhysics , 4077 (1996).[48] V. Kapil, E. Engel, M. Rossi, and M. Ceriotti, Journalof Chemical Theory and Computation , 5845 (2019).[49] B. Cheng, E. A. Engel, J. Behler, C. Dellago, andM. Ceriotti, Proceedings of the National Academy ofSciences , 1110 (2019).[50] Y. Zuo, C. Chen, X. Li, Z. Deng, Y. Chen, J. Behler,G. Csányi, A. V. Shapeev, A. P. Thompson, M. A.Wood, and S. P. Ong, The Journal of Physical Chem-istry A , 731 (2020).[51] A. M. Reilly and A. Tkatchenko, J. Chem. Phys. ,024705 (2013).[52] G. J. O. Beran, Chem. Rev. , 5567 (2016).[53] B. Monserrat, J. G. Brandenburg, E. A. Engel, andB. Cheng, Nature Communications , 5757 (2020).[54] T. Morawietz, A. Singraber, C. Dellago, and J. Behler,Proceedings of the National Academy of Sciences ,8368 (2016).[55] R. Jinnouchi, F. Karsai, and G. Kresse, Physical Re-view B , 014105 (2019).[56] A. Grisafi and M. Ceriotti, The Journal of ChemicalPhysics , 204105 (2019). [57] S. N. Pozdnyakov, M. J. Willatt, A. P. Bartók, C. Or-tner, G. Csányi, and M. Ceriotti, Physical ReviewLetters , 166001 (2020).[58] A. Katrusiak, M. Podsiadło, and A. Budzianowski,Crystal Growth & Design , 3461 (2010).[59] A. Dawson, D. R. Allan, S. A. Belmonte, S. J. Clark,W. I. F. David, P. A. McGregor, S. Parsons, C. R.Pulham, and L. Sawyer, Crystal Growth & Design ,1415 (2005).[60] J.-L. Leviel, G. Auvert, and J.-M. Savariault, ActaCrystallographica Section B: Structural Crystallogra-phy and Crystal Chemistry , 2185 (1981).[61] K. H. Yurtseven and M. Senol, Acta Physica PolonicaA , 698 (2013).[62] P. Lucaioli, E. Nauha, I. Gimondi, L. S. Price, R. Guo,L. Iuzzolino, I. Singh, M. Salvalaglio, S. L. Price, andN. Blagden, CrystEngComm , 3971 (2018).[63] J. G. Brandenburg and S. Grimme, Acta Crystallo-graphica Section B: Structural Science, Crystal Engi-neering and Materials , 502 (2016).[64] J. v. d. Streek and M. A. Neumann, CrystEngComm , 7135 (2011).[65] J. Hoja, A. M. Reilly, and A. Tkatchenko, WIREsComputational Molecular Science , e1294 (2017).[66] J. H. Lloyd-Williams and B. Monserrat, Phys. Rev. B , 184301 (2015).[67] A. Dawson, D. R. Allan, S. A. Belmonte, S. J. Clark,W. I. F. David, P. A. McGregor, S. Parsons, C. R.Pulham, and L. Sawyer, Crystal Growth and Design , 1415 (2005).[68] K. Yamauchi, Shigeki, and K. I. Ando, Journal ofMolecular Structure , 9 (2002).[69] N. F. Francia, L. S. Price, J. Nyman, S. L. Price,and M. Salvalaglio, Crystal Growth & Design , 6847(2020).[70] C. Červinka and G. J. O. Beran, Chemical Science ,4622 (2018).[71] J. Klimeš, The Journal of Chemical Physics ,094506 (2016).[72] T. Gruber, K. Liao, T. Tsatsoulis, F. Hummel, andA. Grüneis, Physical Review X , 021043 (2018).[73] A. Zen, J. G. Brandenburg, J. Klimeš, A. Tkatchenko,D. Alfè, and A. Michaelides, Proceedings of the Na-tional Academy of Sciences , 1724 (2018).[74] R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A.von Lilienfeld, Journal of Chemical Theory and Com-putation , 2087 (2015).[75] C. Schran, K. Brezina, and O. Marsalek, The Journalof Chemical Physics , 104105 (2020).[76] F. Musil, M. J. Willatt, M. A. Langovoy, and M. Ce-riotti, J. Chem. Theory Comput. , 906 (2019).[77] J. Nigam, S. Pozdnyakov, and M. Ceriotti, The Jour-nal of Chemical Physics , 121101 (2020).[78] S. Wengert, G. Csányi, K. Reuter, and J. T. Margraf,Chemical Science (2021), 10.1039/D0SC05765G.[79] J. Behler and M. Parrinello, Phys. Rev. Lett. ,146401 (2007). [80] A. Singraber, J. Behler, and C. Dellago, Journal ofChemical Theory and Computation , 1827 (2019).[81] G. Imbalzano, A. Anelli, D. Giofré, S. Klees, J. Behler,and M. Ceriotti, Journal of Chemical Physics ,241730 (2018).[82] R. K. Cersonsky, B. A. Helfrecht, E. A. Engel, andM. Ceriotti, (2021), arXiv:2012.12253.[83] V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Lit-man, T. Spura, B. Cheng, A. Cuzzocrea, R. H.Meißner, D. M.Wilkins, B. A. Helfrecht, P. Juda, S. P.Bienvenue, W. Fang, J. Kessler, I. Poltavsky, S. Van-denbrande, J. Wieme, and M. Ceriotti, ComputerPhysics Communications , 214 (2018).[84] “DFTB+, a software package for efficient approximatedensity functional theory based atomistic simulations:The Journal of Chemical Physics: Vol 152, No 12,” .[85] T. Krüger, M. Elstner, P. Schiffels, and T. Frauen-heim, The Journal of Chemical Physics , 114110(2005).[86] Y. Eldar, M. Lindenbaum, M. Porat, and Y. Y.Zeevi, IEEE Transactions on Image Processing , 1305(1997).[87] M. Ceriotti, G. A. Tribello, and M. Parrinello, Journalof Chemical Theory and Computation , 1521 (2013).[88] R. J. G. B. Campello, D. Moulavi, A. Zimek, andJ. Sander, ACM Trans. Knowl. Discov. Data , 5(2015).[89] J. P. Perdew, K. Burke, and M. Ernzerhof, PhysicalReview Letters , 3865 (1996).[90] A. Tkatchenko and M. Scheffler, Physical Review Let-ters , 073005 (2009).[91] J. P. Perdew, M. Ernzerhof, and K. Burke, Journal ofChemical Physics , 9982 (1996).[92] C. Adamo and V. Barone, Journal of Chemical Physics , 6158 (1999).[93] A. Tkatchenko, R. A. Di Stasio Jr, R. Car, andM. Scheffler, Physical Review Letters , 236402(2012).[94] A. Ambrosetti, A. M. Reilly, R. A. Di Stasio Jr,and A. Tkatchenko, Journal of Chemical Physics ,018A508 (2014).[95] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu,X. Ren, K. Reuter, and M. Scheffler, ComputerPhysics Communications , 2175 (2009).[96] X. Ren, P. Rinke, V. Blum, J. Wieferink,A. Tkatchenko, A. Sanfilippo, K. Reuter, andM. Scheffler, New Journal of Physics , 053020(2012).[97] S. Levchenko, X. Ren, J. Wieferink, R. Johanni,P. Rinke, V. Blum, and M. Scheffler, ComputerPhysics Communications , 60 (2015).[98] H. J. Monkhorst and J. D. Pack, Physical Review B , 5188 (1976).[99] M. Schlipf and F. Gygi, Computer Physics Communi-cations , 36 (2015).[100] P. Raiteri, J. D. Gale, and G. Bussi, Journal ofPhysics: Condensed Matter , 334213 (2011).[101] D. Frenkel and A. J. C. Ladd, The Journal of ChemicalPhysics81