Communication pubs.acs.org/crystal
Accidental Degeneracy in Crystalline Aspirin: New Insights from High-Level ab Initio Calculations Shuhao Wen and Gregory J. O. Beran* Department of Chemistry, University of California, Riverside, California 92521, United States S Supporting Information *
ABSTRACT: We perform the first high-level ab initio calculations (MP2) on crystalline aspirin using a newly developed fragment-based QM/MM method. Contrary to earlier density functional theory predictions, the two polymorphs are virtually degenerate, which is consistent with experimentally observed intergrowth structures. This near-degeneracy arises “accidentally” from a competition between intramolecular relaxation (form I) and intermolecular hydrogen bonding (form II).
C
interesting molecular crystals like aspirin.19,26 Our hybrid manybody interaction (HMBI) fragment method15−19 combines a QM treatment of individual molecules in the unit cell and their short-range pairwise interactions with a polarizable MM model for long-range and many-body intermolecular interactions:
hanges in crystal packing can significantly alter the bioavailability or other key properties of a pharmaceutical.1 Forty years after the original crystal structure of aspirin (acetylsalicylic acid) was solved,2 a new polymorph of aspirin, form II, was predicted in 20043 and experimentally observed the following year.4 This claim proved highly controversial, however, because the crystal packing differs only slightly between the two polymorphs, and domains of both forms are often found in intergrowth structures.5,6 Only recently has clear evidence for form II as a distinct polymorph emerged,7−9 though a proper explanation of the differential crystal packing interactions remains elusive. The existence of intergrowth structures led to the suggestion that the two forms might be “accidentally” degenerate.6 However, density functional theory (DFT) calculations predict that form II is 2−2.5 kJ/mol more stable than form I,10,11 which is significantly larger than the subkJ/mol energies that separate different packing motifs in other disordered crystals.12−14 Using a new fragment-based hybrid quantum classical method we have developed,15−19 we performed the first high-level electronic structure calculations on crystalline aspirin. We demonstrate that the two aspirin polymorphs are in fact virtually isoenergetic and that this neardegeneracy arises from a competition between intramolecular geometry relaxation and enhanced intermolecular interactions in the crystal. Obtaining reliable theoretical predictions in systems such as aspirin, a flexible molecule whose unit cell contains four symmetry-equivalent molecules and 84 total atoms in the unit cell, is difficult. The model must carefully balance intra- and intermolecular interactions, and conventional electronic structure methods are either too computationally expensive or too unreliable. DFT, for instance, requires corrections for van der Waals dispersion interactions,20−22 and the variations among predictions from different functionals or dispersion corrections can substantially exceed the energy gaps in question.23−25 Newly developed fragment methods now make accurate, high-level electronic structure calculations feasible in chemically © 2012 American Chemical Society
QM MM MM EHMBI = E1QM ‐body + E 2‐body short + E 2‐body long + Emany‐body
(1)
For the MM portions, we utilize an ab initio force field,16,18 which contains terms describing long-range two-body electrostatics, long-range two-body dispersion, self-consistent manybody induction, and three-body Axilrod−Teller dispersion. The force field parameters (distributed multipole moments, distributed polarizabilities, and dispersion coefficients27,28) are computed on-the-fly from quantum mechanical calculations on the individual molecules in the unit cell. Benchmark calculations demonstrate that the HMBI model can predict molecular crystal lattice energies for small molecules to within 1−2 kJ/mol of experiment, or within experimental errors, when very large basis sets and high-level coupled cluster methods are used.18 Error cancellation between two polymorphs can lead to even more accurate predictions for the relative energies, particularly when the packing motifs are similar. In aspirin both polymorphs exhibit stacked layers of hydrogen-bonded carboxylic acid dimers. As shown in Figure 1, the primary difference between forms I and II lies in whether the interlayer hydrogen bonds form dimers (form I) or catemers (form II). To investigate the polymorphism in aspirin, we first optimized the experimental form I and II structures.5,6,29 In principle, we could use HMBI for this purpose, but that would be cost-prohibitive here. DFT structures are usually sufficiently accurate,19,30 so we use the dispersion-corrected B3LYP-D* functional 30 and the TZP basis, as implemented in CRYSTAL09,31,32 to relax the experimental structures. Received: March 17, 2012 Published: April 3, 2012 2169
dx.doi.org/10.1021/cg300358n | Cryst. Growth Des. 2012, 12, 2169−2172
Crystal Growth & Design
Communication
Table 2. Predicted Lattice Energies of Aspirin, in kJ/mol RI-MP2
form I form II ΔEI→II ΔHI→II(0 K)a a
SCS(MI) RI-MP2
aug-ccpVDZ
aug-ccpVTZ
aug-ccpVDZ
aug-ccpVTZ
113.7 113.5 0.2 −0.2
132.1 132.0 0.1 −0.3
132.5 132.3 0.1 −0.3
135.6 135.5 0.0 −0.4
ΔHI→II(0 K) = ΔEI→II + ΔZPEI→II, and ΔZPEI→II = −0.4 kJ/mol.
errors between polymorphs with similar packing leads to consistent predictions of the relative energies that are very insensitive to the user-defined options in the HMBI model (see Supporting Information). Overall, form I is 0.0−0.1 kJ/mol more stable than form II, which is consistent with the accidental degeneracy hypothesis for explaining the experimentally observed polymorph intergrowths. The zero-point energy (ZPE) can play an important role in distinguishing closely spaced polymorphs.24 B3LYP-D* harmonic frequency calculations at the Γ point43,44 indicate that zero-point energy further reverses the nominal polymorph ordering by stabilizing form II 0.4 kJ/mol relative to form I. Still, the two structures remain virtually degenerate at 0 K. Table 3 compares these predictions with those from earlier studies. The MM predictions actually agree best with these
Figure 1. Polymorphs of aspirin. The key difference lies in whether the interlayer hydrogen bonds form dimers (form I) or catemers (form II).
The optimized aspirin structures agree very well with the experimental ones, as shown in Table 1. The calculated lattice parameters and volume lie within 2% of the experimental values, and overlaying the DFT and experimental structures gives rmsd15 values33 of only ∼0.1 Å. These results are consistent with the typical performance of DFT for modeling crystal geometries30 and with the thermal expansion that would occur upon heating from 0 K (DFT) to finite temperatures.35 Additional details are provided as Supporting Information. The most notable difference between the DFT and experimental structures is that form I intramolecular torsion angles for the acetyl group relative to the ring are about 3° smaller than the experimental values. This rotation couples with the compression along the c-axis to decrease the interlayer dimer hydrogen bond distance from 2.6 Å measured at 123 K to 2.4 Å optimized at 0 K. The TZP form II structure exhibits a similar compression of the interlayer hydrogen bonds relative to the 180 K experimental structure, though the catemeric bonding reduces the flexibility of the acetyl side chains. This acetyl flexibility plays an important role in the relative polymorph energetics, as discussed below. Single-point lattice energies were then calculated with HMBI on the DFT geometries (Table 2). The QM calculations were performed using counterpoise-corrected dual-basis resolution of the identity second-order Møller−Plesset perturbation theory (RI-MP2)36,37 and spin-component-scaled MP2 for molecular interactions, SCS(MI) RI-MP2,38 as implemented in Q-Chem 3.1.39 The latter provides a low-cost means of correcting the MP2 overestimation of π-stacking dispersion interactions, which are important in aspirin.17,18,40 The force field parameters were calculated from coupled Kohn−Sham theory with the PBE0 functional41 and the Sadlej basis using Camcasp.42 The predicted lattice energies vary somewhat among the different methods and basis sets, but the large cancellation of
Table 3. Summary of Theoretical Predictions for the Polymorph Energy Differencea method
ΔEI→II (kJ/mol)
source
MM + QM intramolecular COMPASS MM45 B3LYP + dispersion/6-31G** B3LYP-D*/6-31G* B3LYP + dispersion/6-311G** B3LYP-D*/TZP HMBI RI-MP2/aug-cc-pVTZ HMBI SCS(MI) RI-MP2/aug-cc-pVTZ
−0.2 −1.2 −2.5 −2.5 −1.9 −1.4 +0.1 +0.0
ref 3 ref 8 ref 11 (this work) ref 11 (this work) (this work) (this work)
a
All QM results are counterpoise corrected, and zero-point energy is not included.
MP2-level HMBI calculations, but this is largely fortuitous. Pure MM treatments such as COMPASS are known to struggle when trying to predict crystal structures of flexible molecules due to an uneven balance between inter- and intramolecular forces. The original work in ref 3 that first predicted form II compensates for this weakness by using QM-determined intramolecular conformers in their rigid-molecule crystal packing search. Given the small energy differences, however, those authors stated that “little confidence can be placed in the
Table 1. Optimized Lattice Parameters for Aspirin Forms I and II at the B3LYP-D*/TZP Level with Percent Errors Relative to Experiment Listed in Parentheses form I form II
method
a (Å)
b (Å)
c (Å)
β (deg)
vol (Å3)
rmsd15a (Å)
DFT exptb DFT exptc
11.259 (−0.2%) 11.278 12.226 (0.6%) 12.152
6.556 (0.1%) 6.552 6.374 (−2.0%) 6.506
11.080 (−1.7%) 11.274 11.257 (−1.0%) 11.368
96.07 (0.2%) 95.84 111.25 (−0.3%) 111.57
813.3 (−1.9%) 828.7 817.7 (−2.2%) 835.8
0.096 0.114
a
Root-mean-square deviation relative to experiment (excluding hydrogen atoms) from overlaying a 15-molecule coordination sphere of molecules,33 as calculated by Mercury 2.4.34 bRef 5. cRef 6. 2170
dx.doi.org/10.1021/cg300358n | Cryst. Growth Des. 2012, 12, 2169−2172
Crystal Growth & Design
■
relative orderings...”.3 In addition, neither force-field model included polarization, which is necessary to understand the polymorphism here. So while MM methods have proved very useful in many molecular crystal studies, the nominal agreement between MM and HMBI here is at least partially unphysical. As noted earlier, DFT with the B3LYP functional plus a dispersion correction overstabilizes form II by up to 2.5 kJ/mol. The results improve moderately in larger basis sets, but the energy gap is still too large. Overall, the intergrowth structures are indicative of very small energy gaps, and B3LYP-D* appears unable to predict that correctly. To understand this near-degeneracy chemically, we decompose the HMBI energies into their fragment contributions (Figure 2). The intramolecular interactions (1-body terms)
Communication
ASSOCIATED CONTENT
S Supporting Information *
Optimized crystal structures and their comparison with experiment, detailed energetics, and analysis of the uncertainty in the force field contributions. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS Funding for this work from the National Science Foundation (CHE-1112568) and supercomputer time from the Teragrid (TG-CHE110064) are gratefully acknowledged.
■
REFERENCES
(1) Bauer, J.; Spanton, S.; Quick, R.; Quick, J.; Dziki, W.; Porter, W.; Morris, J. Pharm. Res. 2001, 18, 859−866. (2) Wheatley, P. J. J. Chem. Soc. 1964, 6036−6048. (3) Ouvrard, C.; Price, S. L. Cryst. Growth Des. 2004, 4, 1119−1127. (4) Vishweshwar, P.; McMahon, J. A.; Oliveira, M.; Peterson, M. L.; Zaworotko, M. J. J. Am. Chem. Soc. 2005, 127, 16802−3. (5) Bond, A. D.; Boese, R.; Desiraju, G. R. Angew. Chem., Int. Ed. 2007, 46, 615−7. (6) Bond, A. D.; Boese, R.; Desiraju, G. R. Angew. Chem., Int. Ed. 2007, 46, 618−22. (7) Chan, E. J.; Welberry, T. R.; Heerdege, A. P.; Goossens, D. J. Acta Crystallogr., B 2010, 66, 696−707. (8) Bauer, J. D.; Haussuhl, E.; Winkler, B.; Arbeck, D.; Milman, V.; Robertson, S. Cryst. Growth Des. 2010, DOI: 100527101541018. (9) Bond, A. D.; Solanko, K. A.; Parsons, S.; Redder, S.; Boese, R. CrystEngComm 2011, 13, 399. (10) Li, T.; Feng, S. Pharm. Res. 2006, 23, 2326−2332. (11) Li, T. J. Pharm. Sci. 2007, 96, 755−760. (12) Copley, R. C. B.; Barnett, S. A.; Karamertzanis, P. G.; Harris, K. D. M.; Kariuki, B. M.; Xu, M.; Nickels, E. A.; Lancaster, R. W.; Price, S. L. Cryst. Growth Des. 2008, 8, 3474−3481. (13) Torrisi, A.; Leech, C. K.; Shankland, K.; David, W. I. F.; Ibberson, R. M.; Benet-Buchholz, J.; Boese, R.; Leslie, M.; Catlow, C. R. a.; Price, S. L. J. Phys. Chem. B 2008, 112, 3746−58. (14) Winkel, K.; Hage, W.; Loerting, T.; Price, S. L.; Mayer, E. J. Am. Chem. Soc. 2007, 129, 13863−71. (15) Beran, G. J. O. J. Chem. Phys. 2009, 130, 164115. (16) Sebetci, A.; Beran, G. J. O. J. Chem. Theory Comput. 2010, 6, 155−167. (17) Beran, G. J. O.; Nanda, K. J. Phys. Chem. Lett. 2010, 1, 3480− 3487. (18) Wen, S.; Beran, G. J. O. J. Chem. Theory Comput. 2011, 7, 3733−3742. (19) Wen, S.; Nanda, K.; Huang, Y.; Beran, G. J. O. Phys. Chem. Chem. Phys. 2012, DOI: 10.1039/c2cp23949c. (20) Grimme, S. WIRES: Comput. Mol. Sci. 2011, 1, 211−228. (21) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (22) Lu, D.; Li, Y.; Rocca, D.; Galli, G. Phys. Rev. Lett. 2009, 102, 206411. (23) Hongo, K.; Watson, M. A.; Sanchez-Carrera, R. S.; Iitaka, T.; Aspuru-Guzik, A. J. Phys. Chem. Lett. 2010, 1, 1789−1794. (24) Rivera, S. A.; Allis, D. G.; Hudson, B. S. Cryst. Growth. Des. 2008, 8, 3905−3907. (25) Karamertzanis, P. G.; Day, G. M.; Welch, G. W. A.; Kendrick, J.; Leusen, F. J. J.; Neumann, M. A.; Price, S. L. J. Chem. Phys. 2008, 128, 244708.
Figure 2. Energy decomposition of the SCS(MI) RI-MP2/aug-ccpVTZ polymorph energy difference into the QM, MM (2-body electrostatics, induction, 2-body dispersion, and 3-body dispersion), and ZPE contributions.
favor form I by 1.5 kJ/mol, most of which comes from differences in the rotation of the acetyl group. The interlayer dimer structure of form I allows the flexible acetyl groups to relax somewhat, while the catemeric structure of form II constrains the side-chain in a higher-energy conformer. On the other hand, form II exhibits more attractive intermolecular interactions. Most of the additional intermolecular stabilization of form II comes from the longer-range two-body electrostatics and the many-body induction contributions. The catemer chains between layers found in form II stabilize it relative to form I through hydrogen-bond cooperativity. In summary, we have demonstrated that forms I and II of aspirin are virtually isoenergetic, which explains why intergrowths of the two packing motifs frequently occur. This “accidental” degeneracy arises from a competition between a more favorable intramolecular conformation (form I) with improved intermolecular hydrogen-bond cooperativity effects (form II). This study also exemplifies how recent advances in fragment-based electronic structure methods create new opportunities for obtaining insights into chemically interesting molecular crystals. 2171
dx.doi.org/10.1021/cg300358n | Cryst. Growth Des. 2012, 12, 2169−2172
Crystal Growth & Design
Communication
(26) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. Chem. Rev. 2011, 112, 632−672. (27) Stone, A. J. The Theory of Intermolecular Forces; Clarendon Press: Oxford, 2002; pp 5−10. (28) Stone, A. J.; Misquitta, A. J. Int. Rev. Phys. Chem. 2007, 26, 193− 222. (29) Cambridge Crystallographic Data Center reference codes ACSALA14 (form I) and ACSALA15 (form II). (30) Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. CrystEngComm 2008, 10, 405−410. (31) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. Z. Kristallogr. 2005, 220, 571−573. (32) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M.; Science, C.; Technologies, A. CRYSTAL09 User’s Manual; University of Torino: Torino, 2009. (33) Chisholm, J. A.; Motherwell, W. D. S. J. Appl. Crystallogr. 2005, 38, 228−231. (34) Macrae, C. F.; Bruno, I. J.; Chisholm, J. A.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Rodriguez-Monge, L.; Taylor, R.; van de Streek, J.; Wood, P. A. J. Appl. Crystallogr. 2008, 41, 455−470. (35) Beyer, T.; Price, S. L. CrystEngComm 2000, 2, 183. (36) Steele, R. P.; Distasio, R. A.; Shao, Y.; Kong, J.; Head-Gordon, M. J. Chem. Phys. 2006, 125, 074108. (37) Steele, R. P.; Distasio, R. A.; Head-Gordon, M. J. Chem. Theory Comput.. 2009, 5, 1560−1572. (38) Distasio, R. A.; Head-Gordon, M. Mol. Phys. 2007, 105, 1073− 1083. (39) Shao, Y.; et al. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191. (40) Riley, K. E.; Pitonák, M.; Jurecka, P.; Hobza, P. Chem. Rev. 2010, 110, 5023−63. (41) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (42) Misquitta, A. J.; Stone, A. J. CamCASP v5.6 (2011), http:// www-stone.ch.cam.ac.uk/programs.html. Accessed Februrary 23, 2011. (43) Pascale, F.; Zicovich-Wilson, C. M.; López Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 888−97. (44) Zicovich-Wilson, C. M.; Pascale, F.; Roetti, C.; Saunders, V. R.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 1873−81. (45) Sun, H. J. Phys. Chem. B 1998, 102, 7338−7364.
2172
dx.doi.org/10.1021/cg300358n | Cryst. Growth Des. 2012, 12, 2169−2172
