In this paper, we report a novel crystal structure prediction method and demonstrate its accuracy on a large set of diverse molecules. The method integrates a novel systematic approach to search the crystal packing parameters and a hierarchical energy ranking method that balances accuracy and cost (Figure ). The new packing search method uses a divide-and-conquer strategy to break down the parameter space into subspaces based on space group symmetries. Each subspace is then searched consecutively. The energy ranking method combines molecular dynamics (MD) simulations using a classical force field (FF), structure optimization and reranking using a machine learning force field (MLFF) with long range electrostatic and dispersion interactions, and periodic density functional theory (DFT) calculations for ranking the final shortlist. The temperaturedependent stability of different polymorphs is evaluated with free energy calculations using previously established methods. Currently focused on searching crystal structures with one molecule in the asymmetric unit (ASU), corresponding to the Z'=1 search space, the method was validated on a large set of 66 molecules with 137 unique crystal structures, including all relevant molecules from the first six CCDC CSP blind tests, Target XXXI from the seventh blind test, other molecules studied by previous CSP methods, and several molecules from modern drug discovery programs. For all the molecules in this large test set, the experimentally known polymorphs are correctly predicted by our method and are ranked among the top candidate structures. For several molecules, our prediction suggests new low energy polymorphs yet to be discovered by experiment, implying potential risks that could jeopardize the development of the currently known forms of these compounds. Comparisons to the results of other computational approaches are made when feasible.