Figure Overview of the computational polymorph prediction method. The novel packing search method uses a divide-and-conquer strategy to break down the parameter space into subspaces based on space group symmetries. 3D candidate crystal structures are built step by step from low energy conformers to different sized clusters with favorable interactions following the funnel shaped energy landscape mimicking the crystal nucleation process. Each conformer can generate multiple candidate structures via different pathways (the entire sampling tree), and multiple pathways can lead to the same candidate structure (the blue colored paths generate the same candidate structure). The energy ranking method integrates a multi-stage energy relaxation and filtering process with increasing accuracy and cost, including molecular dynamics simulations using the classical force field, structure optimization and reranking using a machine learning force field (MLFF) with long range electrostatic and dispersion interactions, and periodic DFT calculations for ranking the final shortlist. For cocaine (a), the only known experimental structure is predicted to have much lower energy than all other predicted candidate structures with an energy gap of about 1.0 kcal/mol by PBD-D3 and PBE0-MBD functionals (e). For MK-2022 all functionals predict experimentally known form I to be most stable. For olanzapine, the three resolved experimental structures are predicted among the top five candidate structures with three different DFT functionals (f). For the unresolved form III, two predicted candidate structures have the spectra characteristics of form III from experiment, one corresponding to a previously identified structure in the Pbca space group (j), and one novel structure in the P21/c space group with an identical 2D layer (i). For GSK269984B, candidate structures with favorable intermolecular hydrogen bonds (k) and (l), including the typical bidentate hydrogen bonds for carboxylic acids (k), are predicted to have lower energy than the experimental structure with only intramolecular hydrogen bonds (m), indicating the possibility of other polymorphs with competitive stability yet to be discovered by experiment. ). For (a) ROY and (b) Galunisertib, molecules holding the current record for the largest number of fully characterized polymorphs, our CSP approach correctly predicted all the known Z'=1 polymorphs. For rotigotine (c), where the late appearance of the more stable form II disrupted the clinical usage of the original form I formulation, our CSP calculations correctly predicted all the experimentally known forms I and II, the candidate form III from a prior CSP study, and another novel candidate crystal in the P212121 space group with greater stability than form I. For MK-8876 (d) and MK-2202 (e), a previous CSP study using the GRACE software correctly predicted a candidate structure matching the experimentally known polymorph for MK-2022, but failed for MK-8876. Our CSP calculations correctly predicted both experimental structures with the experimental structure ranking 1st for MK-2022 and 2nd for MK-8876. For LY156735 (f), a potent and selective melatonin agonist with some unusual crystallization behavior among the R and S enantiomers, our calculations correctly predicted all three known experimental forms, and suggested a few other competitive structures, indicating the complex polymorphic landscape of this molecule may remain to be extensively explored experimentally. The predicted lowest energy structure matching the experimental Z'=1 structure of Mod I is most stable at all temperatures below 400K according to the calculations, matching experimental results that Mod I is the most stable form at all temperatures below its melting point of 361 K.