Surface defects can significantly affect the performances of energetic crystals, particularly safety and mechanics, and deepening the insight of the surface defect formation mechanism is expected to facilitate the defect engineering. Nevertheless, the mechanism generally remains unclear, primarily due to a long time span of the formation, making it difficult to capture it by conventional molecular simulations. To address this, this work takes the most powerful energetic compound commercialized already, 2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane (CL-20), as a sample; combines quantum chemistry and machine learning methods to refine a force field for it, the OPLS-CL-20 force field; and performs well-tempered metadynamics simulations on the main crystal faces of ε-CL-20 in six solvents to explore the vacancy defect formation mechanism thereby. It is found that, among all the main crystal faces, (002) exhibits the largest molecular detachment work across all solvents and preserves the highest perfection, and all faces show greater molecular detachment work in the DMF solution relative to others, implying the highest possibility of ε-CL-20 keeping perfect therein. Our study may provide theoretical guidance for the solvent selection and crystal quality control of molecular crystalline materials.