In this study, we present a systematic computational investigation to analyze the long-debated free energy stability of two wellknown aspirin polymorphs, denoted as Form I and Form II. Specifically, we developed a strategy to collect training configurations covering diverse interatomic interactions between representative functional groups in aspirin crystals. Utilizing a state-of-the-art neural network interatomic potential (NNIP) model, we trained an accurate machine learning potential to simulate aspirin crystal dynamics under finite temperature conditions with ∼0.46 kJ/mol/molecule accuracy. Employing the trained NNIP model, we performed thermodynamic integration to assess the free energy difference between aspirins I and II, accounting for the anharmonic effects in a large supercell consisting of 512 molecules. For the first time, our results convincingly demonstrated that Form I is more stable than Form II at 300 K, ranging from 0.74 to 1.83 kJ/mol/molecule, aligning with experimental observations. Unlike the majority of previous simulations based on (quasi)harmonic approximations in a small super cell, which often found degenerate energies between aspirins I and II, our findings underscore the importance of anharmonic effects in determining polymorphic stability ranking. Furthermore, we proposed the use of the rotational degrees of freedom of methyl and ester/phenyl groups in aspirin crystals as characteristic motions to highlight rotational entropic contribution that favors the stability of Form I. From the structural perspective, we also found that the subtle free energy difference can be used to explain the distinct thermal expansion responses as observed in both experimental and simulation data. Beyond the aspirin polymorphism, we anticipate that such entropy-driven stabilization can be broadly applicable to many other organic systems, suggesting that our approach holds great promise for stability studies in small-molecule drug design.
