Crystal polymorphism of complex liquids plays a crucial role in industrial crystallization, food technology, pharmaceuticals, and materials engineering. However, the experimental identification of unknown crystal structures can be challenging, particularly for high-viscosity complex liquids, such as ionic liquids (ILs). In this study, we performed a molecular dynamics simulation coupled with metadynamics to investigate an imidazolium IL (1-alkyl-3-methylimidazolium hexafluorophosphates). The simulation employed two distinct radial-distribution functions, represented by Gaussian window functions as collective variables, and revealed at least two crystal-like phases distinct from the known α and β crystal phases typically formed by this IL. Additionally, the simulation unveiled a unique phase characterized by the ordered spatial arrangement of anion aggregations. These crystal-like and unique phases emerged regardless of the potential used. The simulation methodology presented here is broadly applicable for exploring unknown phases in complex systems and contributes to the design of functional materials, such as porous ILs for gas molecule capture and separation.