  In this study, we present a systematic computational investigation to analyze
the long debated crystal stability of two well known aspirin polymorphs,
labeled as Form I and Form II. Specifically, we developed a strategy to collect
training configurations covering diverse interatomic interactions between
representative functional groups in the aspirin crystals. Utilizing a
state-of-the-art neural network interatomic potential (NNIP) model, we
developed an accurate machine learning potential to simulate aspirin crystal
dynamics under finite temperature conditions with $\sim$0.46 kJ/mol/molecule
accuracy. Employing the trained NNIP model, we performed thermodynamic
integration to assess the free energy difference between aspirin Forms 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 the experimental observations. Unlike the
majority of previous simulations based on (quasi)harmonic approximations in a
small super cell, which often found the degenerate energies between aspirin I
and II, our findings underscore the importance of anharmonic effects in
determining polymorphic stability ranking. Furthermore, we proposed the use of
rotational degrees of freedom of methyl and ester/phenyl groups in the aspirin
crystal, as characteristic motions to highlight rotational entropic
contribution that favors the stability of Form I. Beyond the aspirin
polymorphism, we anticipate that such entropy-driven stabilization can be
broadly applicable to many other organic systems and thus our approach,
suggesting our approach holds a great promise for stability studies in small
molecule drug design.
