Our CSP workflow begins by converting a molecule's bonding information to a single, canonicalized three dimensional geometry. This geometry is input to the Schrödinger Force Field Builder (FFBuilder), which examines the coverage of the OPLS4 force-field for the molecule and refits torsion parameters which are not well represented in the default training set with additional QM calculations when needed. This approach differs from other CSP approaches where customized parameters are trained specifically on one molecule or a set of molecules. The resulting parameters from our method are transferable to different chemistries and for different applications, whereas the customized parameters used in other CSP methods only have limited scope of application. More detailed discussion on the construction of OPLS4 and automated torsion parameter fitting is reported elsewhere. To generate a pool of conformers for molecular crystal packing, the molecule is fragmented, low energy ring states are determined, pyramidal nitrogens with unique inversion states are identified, and energy profiles of each flexible torsion are obtained. Sample angles for each flexible torsion are determined using a combination of simple chemical rules and automated detection of low energy regions. All combinations of the chosen torsion angles, ring states, and pyramidal nitrogen states are placed onto the molecule fragments. Fragments with low energy according to the OPLS4 force field are combined to produce a diverse pool of unrelaxed conformers. The conformers within 8 kcal/mol of the global minimum are then relaxed, combined with low-energy unrelaxed conformers, and deduplicated to produce the conformer pool. Our crystal packing search treats the sampled conformers as rigid building blocks and only targets crystals with Z'=1 currently. Molecular clusters constructed from initial guesses of unevenly distributed parameters corresponding to basic symmetry elements, such as translation, inversion, rotation, screw, and glide are minimized against these parameters. Additional symmetry elements are searched sequentially for each cluster until unit cells can be constructed for the space group under consideration. The cluster size ranges from 2 for an inversion dimer to 36 for the final representation of the 3D crystal in some space groups. For example, the P1 space group contains three translation symmetries, and we effectively search low-energy periodic structures in one, two, and three dimensions consecutively corresponding to clusters with 3, 6, and 12 molecules (Figure (a)). When multiple different symmetries exist, the search order matters. Therefore we search with different symmetry orderings for these space groups. The search scope is determined by the energy tolerance of the symmetrized clusters at each search step, and intermediate clusters with energies lower than 6 kcal/mol per molecule from the global minimum of that step are retained for the following steps. As a result, this systematic and hierarchical approach does not need a feedback loop to determine a termination point, unlike other CSP methods based on random sampling. After the packing search, the predicted crystals pass through a multi-stage energy relaxation and filtering of increasing accuracy and cost, including MD, QRNN and DFT stages where all structural degrees of freedom are fully flexible to relax. At each stage, high energy configurations are discarded. The crystal candidates first go through structural relaxation and filtering using symmetry constrained molecular dynamics with the OPLS4 force field. Then, the machine learning force field (MLFF) is used for relaxation and reranking of the remaining crystals. Following this, a coarse DFT relaxation using PBE-D3 is performed followed by more DFT evaluations with increasingly accurate settings and functionals to obtain the final energy ranking. On average, it takes about 20 GPU seconds to relax a crystal with MD, 20 CPU minutes with QRNN, and 60 CPU hours with DFT.