Role of Dispersion Interactions in the Polymorphism and Entropic Stabilization of the Aspirin Crystal

Aspirin has been used and studied for over a century but has only recently been shown to have an additional polymorphic form, known as form II. Since the two observed solid forms of aspirin are degenerate in terms of lattice energy, kinetic effects have been suggested to determine the metastability of the less abundant form II. Here, first-principles calculations provide an alternative explanation based on free-energy differences at room temperature. The explicit consideration of many-body van der Waals interactions in the free energy demonstrates that the stability of the most abundant form of aspirin is due to a subtle coupling between collective electronic fluctuations and quantized lattice vibrations. In addition, a systematic analysis of the elastic properties of the two forms of aspirin rules out mechanical instability of form II as making it metastable. The ability of molecules to yield multiple solid forms, or polymorphs, has significance for diverse applications ranging from drug design and food chemistry to nonlinear optics and hydrogen storage [1]. For example, different solid forms of an active pharmaceutical ingredient can affect its bioavailability and formulation, sometimes in unpredictable ways [2,3]. The computational modeling of polymorphism has seen many advances in recent years, in both generation of reasonable structures [4,5] and optimization and accurate ranking of these polymorphs [6,7]. Even when polymorphs can be predicted for a given molecule, theory and experiment often struggle to understand why a given polymorph is stable under certain thermodynamic conditions [8]. Aspirin (acetylsalicylic acid) is a widely used analgesic that clearly illustrates many of the challenges of understanding polymorphism. A second less common form of aspirin has only been predicted [9] and characterized [10,11] in recent years. State-of-the-art quantum-chemical calculations predict very small energy differences of the order of AE0.1 kJ=mol between the two polymorphs [12,13], making both solid forms essentially degenerate in terms of lattice energy. This raises the question of why the second form took so long to be discovered and why the first form appears to be more abundant. A possible explanation is that kinetic effects, such as slow growth of form II, may play a role, and, indeed, growth of form II can promoted by certain conditions [10,11]. Nanoindentation experiments suggest that, despite their similar structures (Fig. 1), the two polymorphs have markedly different mechanical properties, with form II appearing to be softer and potentially susceptible to shear instability [14]. However, computational studies of their elastic …

The ability of molecules to yield multiple solid forms, or polymorphs, has significance for diverse applications ranging from drug design and food chemistry to nonlinear optics and hydrogen storage [1].For example, different solid forms of an active pharmaceutical ingredient can affect its bioavailability and formulation, sometimes in unpredictable ways [2,3].The computational modeling of polymorphism has seen many advances in recent years, in both generation of reasonable structures [4,5] and optimization and accurate ranking of these polymorphs [6,7].Even when polymorphs can be predicted for a given molecule, theory and experiment often struggle to understand why a given polymorph is stable under certain thermodynamic conditions [8].
Aspirin (acetylsalicylic acid) is a widely used analgesic that clearly illustrates many of the challenges of understanding polymorphism.A second less common form of aspirin has only been predicted [9] and characterized [10,11] in recent years.State-of-the-art quantum-chemical calculations predict very small energy differences of the order of AE0.1 kJ=mol between the two polymorphs [12,13], making both solid forms essentially degenerate in terms of lattice energy.This raises the question of why the second form took so long to be discovered and why the first form appears to be more abundant.A possible explanation is that kinetic effects, such as slow growth of form II, may play a role, and, indeed, growth of form II can promoted by certain conditions [10,11].Nanoindentation experiments suggest that, despite their similar structures (Fig. 1), the two polymorphs have markedly different mechanical properties, with form II appearing to be softer and potentially susceptible to shear instability [14].However, computational studies of their elastic properties have also given conflicting results [9,15].To accentuate this controversy, form II has been observed to revert to form I slowly at room temperature and upon grinding [14], suggesting that form I is thermodynamically more stable.However, no viable mechanism for its thermodynamic stability has yet been established.
Recently, a number of studies have highlighted the importance of many-body van der Waals (vdW) interactions in condensed molecular systems, especially in the context of vdW-inclusive density-functional theory (DFT) [7,[16][17][18].The many-body dispersion (MBD) approach [19,20] has been shown to systematically improve the accuracy of DFT in predicting absolute and relative stabilities of molecular crystals [7,17,20,21].In this contribution, we apply the DFT þ MBD method to aspirin, showing that the role of many-body vdW interactions goes beyond lattice stabilities and energies to affect vibrational and elastic response properties, revealing an entropically driven mechanism behind the relative stability of form-I aspirin.
We begin by putting our approach in the context of other calculations of the relative stabilities of aspirin.These calculations have typically considered the energy differences between the computationally optimized, 0 K lattice of the two polymorphs.Table I gives the MBD lattice energy difference, when coupled with the nonempirical Perdew-Burke-Ernzerhof (PBE) density functional [22] (see the Supplemental Material [23] for further details).For comparison, the value with the pairwise Tkatchenko-Scheffler (TS) scheme for vdW interactions [28] is also reported.Both the MBD and TS results confirm the picture of nearly degenerate polymorphs in terms of lattice energy, already firmly established by quantum-chemical calculations [12].The same conclusion holds when the more accurate hybrid PBE0 functional is used [29].PBE0 þ MBD calculations give very similar results to PBE þ MBD, differing by only 0.1 kJ=mol.
While the two aspirin structures are degenerate in terms of lattice energies, at finite temperature free-energy differences dictate the thermodynamic stability and driving force for the formation of a polymorph.At the moderate temperatures of interest in this work, the free energy of a molecular crystal can be determined from its lattice dynamics in the harmonic approximation [30] FðTÞ where the integral goes over all frequencies ω, gðωÞ is the phonon density of states, E elec is the DFT total energy, and E ZPE is the zero-point energy.Here we report the PBE þ TS and PBE þ MBD free-energy differences at room temperature (298 K) in Table I, along with ZPE corrected totalenergy differences at 0 K.The ZPE contributions are minimal, amounting to 0.2-0.3kJ=mol, in line with previous calculations (0.4 kJ=mol) [12].However, for the free-energy differences at room temperature, the pairwise and many-body descriptions of vdW interactions give qualitatively different pictures.PBE þ TS predicts a slightly more stable form II, with a minor vibrational contribution of 0.5 kJ=mol.In contrast, PBE þ MBD predicts form I to be more stable by the more substantial amount of 2.56 kJ=mol.
We have noted previously that many-body vdW interactions can have a substantial impact on absolute lattice energies [21], including that of form-I aspirin [17].It is demonstrated here that many-body vdW interactions can also substantially affect vibrational properties.The majority of the PBE þ MBD free-energy difference arises from thermal vibrational contributions.To understand this, we plot the phonon density of states of aspirin forms I and II, with the pairwise and many-body description of vdW interactions in Fig. 2. Many-body vdW interactions lead to an appearance of low-frequency phonon modes near 30 cm −1 (four phonon modes per unit cell from 28 to 34 cm −1 ) and additional shifting of phonon frequencies of form I to lower values below 200 cm −1 .Some shifting of the form-II frequencies occurs but to a much lesser extent, FIG. 2 (color online).Phonon density of states for forms I (top) and II (bottom) of aspirin with the DFT-PBE functional in combination with a pairwise (TS) and many-body (MBD) description of vdW interactions.Gaussian broadening of 10 cm −1 is used in this plot.The shaded region highlights the difference in low-frequency phonon modes.TABLE I.The relative total-energy differences (ΔE I→II ¼ E II − E I ) between the two forms of aspirin as determined by using PBE þ TS and PBE þ MBD, together with zero-point energy (ZPE) corrected differences and free-energy differences.
especially in the ≤ 100 cm −1 region.The lower-frequency vibrations of form I are more easily thermally populated, increasing its entropy and ultimately leading to a lower free energy.We stress that the phonon peak at 30 cm −1 in Fig. 2 can be unequivocally attributed to vdW interactions, since this peak is completely absent in PBE phonon calculations using the optimized PBE þ MBD crystal structure (see the Supplemental Material [23]).Experimental terahertz spectra that exhibit low-frequency modes near 37 cm −1 in aspirin pellets at 77 K [31] provide further support for our findings.
Figure 3(a) shows one of the low-frequency vibrations of form I obtained in PBE þ MBD calculations.A number of these modes involve a mixture of motions of whole molecules with concerted motions of methyl groups from different molecules in and out of phase with one another.Form II has comparable modes featuring similar motion of the methyl groups but at higher frequencies [Fig.3(b)].Close inspection of the structure reveals that, while the MBD and TS geometries of the molecule are very similar, the methyl groups are oriented differently with the pairwise and many-body vdW methods.The torsion indicated in Fig. S1 in the Supplemental Material [23] is 47.1°with PBE þ TS.The PBE þ MBD value of 60.1°is in much better agreement with the low-temperature (20 K) neutrondiffraction value of 58.2° [32].A similar difference is seen for the methyl group in form II with the difference in the MBD and TS torsion also being more than 10°.
As the internal structural change is of similar magnitude in both forms, it alone cannot explain the different phonon responses in form I and form II. To explain the lowfrequency phonon peak at 30 cm −1 in form I, we recall that the MBD energy is expressed as [19] where N determines the number of atoms in the system, ω i are the collective dipole excitation frequencies of the whole system (plasmonlike neutral excitations), and ω 0;i are the effective atomic excitation frequencies.We remark that ω i values obtained from MBD along with the corresponding eigenvectors (oscillator strengths) quantitatively reproduce the dielectric constants for a variety of molecular and semiconductor crystals [19,33].By construction, the values of ω 0;i in the MBD method depend only on the atom type and not on the geometry of the system.However, the static polarizability that determines the coupling between the atoms is sensitive to the geometry and leads to renormalization of the ω i frequencies.We note that the characteristic frequencies of the lowest phonon modes in molecular crystals (tens of meV) are 1000 times smaller than typical plasmon energies (8 eV and higher in aspirin).However, the change of the vdW energy associated to characteristic phonon motions [ P 3N p¼1 ðω initial p − ω final p Þ for two aspirin crystal geometries] is on the order of several meV, making plasmon-phonon coupling possible and leading to a visible modification of phonon frequencies in form I, as shown in Fig. 2.
The different hydrogen-bonding motifs in the two forms [14,34] may explain why significantly more softening (and hence plasmon-phonon coupling) is seen for form I. In form I, there are cyclic dimer C − H Á Á Á O hydrogen bonds between the layers of aspirin dimers, leading to an equivalent separation of 4.5 Å (PBE þ MBD geometry) for the methyl C atoms shown in Fig. 3(a) (note that they are not part of the same C − H Á Á Á O dimer), while in form II catemer hydrogen bonding occurs, leading to the methyl groups being separated by 3.74 and 5.51 Å [Fig.3(b)].This asymmetry in methyl-methyl interactions in form II and its strong hydrogen bonds [12] limits the extent of coupling between plasmons and the phonons involving concerted motions of the methyl groups.Comparison of the two phonon modes in the two forms (as in Fig. 3) illustrates how the asymmetry affects the motion of the methyl groups relative to one another.The methyl groups of form I are know to librate substantially at higher temperatures, as has been seen with neutron diffraction [32].Neutron diffraction or inelastic neutron scattering of form II might reveal more insight into the differences between the methyl-group motions and the influence of the C − H Á Á Á O hydrogen bonding on phonons in the two forms.
The free-energy differences between the aspirin polymorphs obtained with PBE þ MBD provide a meaningful explanation for finite-temperature thermodynamic preference for form I. Kinetic factors and mechanical instability of form II have also been suggested as possible reasons for form II appearing to be less stable and common [9][10][11]14].Early computational work suggested that form II might be mechanically unstable [9], which has been supported by nanoindentation experiments on certain faces of the two forms [14] but challenged by semilocal DFT calculations (without any vdW interactions considered) [15].
Here we report the bulk and shear moduli of the two forms calculated with the PBE þ TS and PBE þ MBD  II.We first note that MBD reduces the bulk modulus by over 20% compared to TS, although the experimental value is even smaller, quite likely due to vibrational contributions [35].The two forms have similar bulk moduli, but both the pairwise and many-body descriptions of vdW interactions give form II having a larger modulus, making it less compressible.The level of the vdW treatment affects the qualitative ordering of the shear modulus, with MBD predicting form I to have the lower shear modulus.For a crystal to be stable, the eigenvalues of the elastic-constant matrix must all be positive, satisfying the Born-stability criterion [36].Both forms have all positive eigenvalues with the lowest values being λ min ¼ 2.60 GPa for form I and λ min ¼ 1.86 GPa for form II. At 300 K, λ min ¼ 2.40 GPa experimentally for form I, suggesting that thermal contributions will not substantially affect these values.Thus, both crystals are expected to be stable to compression and shear at room temperature.
Experimental nanoindentation results have suggested that form II is "considerably softer" [14] than form I, supporting the idea of a mechanically unstable form II.However, these experiments are limited to considering the Young's modulus in directions perpendicular to only a few specific faces of the crystals.Figure 4 shows the full spherical plots of the Young's moduli of forms I and II calculated with PBE þ MBD.The Young's modulus is quite complex, illustrating the directionality of the different interactions present in the crystal.Close inspection will show that along certain directions form II is indeed "softer," but that form II has a much larger maximum value for the Young's modulus in the e 1 e 3 plane.The nanoindentation experiments did not probe along this "hard" direction for form II.This direction corresponds to compressing transverse to many of the hydrogen bonds, which are thought to be stronger in form II [12].Comparison between the PBE þ MBD and PBE þ TS Young's moduli for both forms shows qualitative differences, further highlighting the role of many-body vdW interactions in elastic properties.
In the present contribution, we have shown that firstprinciples calculations predict form I of aspirin to be more thermodynamically stable than form II, in agreement with a number of experimental observations of the metastable nature of form II. Crucially, this requires considering many-body van der Waals interactions, not only in lattice energies but also in the lattice vibrations that dictate the finite-temperature vibrational contributions to the free energies.We uncover a novel coupling mechanism between low-frequency phonon modes and long-wavelength collective electron fluctuations, which shifts phonon frequencies in form I to lower values.An analysis of the elastic properties of the two forms has also been performed.Many-body vdW interactions again alter the quantitative and qualitative picture of these results, reducing bulk and shear moduli by more than 20%.Overall, the elastic analysis shows that both forms are expected to be mechanically stable.The more general conclusion of our study is that low-energy vibrational modes of molecular and hybrid  organic-inorganic [37] systems may be significantly modified upon accurate inclusion of vdW interactions in firstprinciples calculations.

FIG. 3 (
FIG. 3 (color online).A low-frequency vibrational mode in(a) form-I aspirin as calculated with PBE þ MBDðω ¼ 33 cm −1 Þ and (b) form-II aspirin as calculated with PBE þ MBDðω ¼ 64 cm −1 Þ.The vectors indicate the direction and relative magnitude of each atom's motion.The black lines indicate the unit cells.

FIG. 4 (
FIG. 4 (color online).Spherical plots of the Young's modulus of form I and form II aspirin as determined by using DFT-calculated elastic constants.For form I, the e 2 and e 3 axes are parallel to the b and c lattice vectors, with e 1 being perpendicular to e 2 and e 3 .To aid comparison, the form-II elastic constants have been rotated around the e 1 axis to make the common, central aspirin dimers of the two forms overlap (cf.Fig.1).

TABLE II .
[15]bulk modulus (B V ) and shear modulus (G V ) of both forms of aspirin determined by DFT calculations.Experimental values for form I were calculated by using elastic constants from Ref.[15].