CRYSTAL GROWTH & DESIGN 2004 VOL. 4, NO. 6 1119-1127
Articles Toward Crystal Structure Prediction for Conformationally Flexible Molecules: The Headaches Illustrated by Aspirin Carole Ouvrard and Sarah L. Price* Department of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom Received February 26, 2004
ABSTRACT: A crystal structure prediction study was carried out on aspirin, based on an analysis of its gas phase conformers and multiple searches for minima in the lattice energy with the molecule held rigid in low energy conformations. Various high levels of ab initio theory were used to estimate the gas phase conformations and energy differences, and accurate distributed multipole-based electrostatic models were used to estimate the electrostatic contribution to the lattice energies. The molecular conformation adopted in the crystal structure is close to a local minimum found in the gas phase ab initio energy using a B3LYP/6-31G(d,p) calculation. A MP2 optimization gives larger differences from the solid state molecular structure. The calculation using the B3LYP molecular conformer predicts the observed crystal structure as one of the most thermodynamically stable and generally the most plausible crystal structure. Alternative molecular conformers, including the gas phase global minimum energy structure and planar transition states, were shown to give less favorable lattice energies. In contrast with a previous study with a flexible molecule force field, the possibility of a planar conformer of aspirin in a crystal structure appears to be unlikely because of the intramolecular energy cost. Although the use of ab initio optimized molecular conformers clearly shows promise for crystal structure prediction for some flexible molecules, the sensitivity of the lattice energies to small distortions of the molecular conformation shows that it can only be used with care when the conformational energy profile is suitable. Introduction Polymorphism,1
the phenomenon that an organic molecule can crystallize in more than one structure, is a major challenge to our fundamental understanding of crystallization and is of immense practical importance in the quality control of the molecular materials industries, such as pharmaceuticals. Conformational polymorphism, where the molecule has very different conformations in its different crystal structures, is quite common. Six different conformations of 5-methyl-2-[(2nitrophenyl)amino]-3-thiophenecarbonitrile are found in different crystal structures,2 giving polymorphs with distinctive colors and morphologies. Abbott Laboratories had to rapidly reformulate the anti-HIV drug Ritonavir3 when a more stable polymorph suddenly appeared in an established manufacturing process, with a different conformation of the molecule producing the more stable form.4 The field of computational crystal structure prediction has tended to concentrate on rigid molecules. The * To whom correspondence should be addressed. E-mail: s.l.price@ ucl.ac.uk.
problems of considering the huge range of possible crystal structures, and modeling the intermolecular forces involved, are already hugely challenging, without the extra degrees of freedom and conformational energies required for flexible molecules. Most of these crystal structure predictions are based on searching for the global minimum in the lattice energy, Ulatt. Despite the crudity of this nevertheless computationally demanding method, there have been some successful crystal structures predicted under blind test conditions.5,6 The limitationsofthisapproachhavebeenwidelydiscussed,7-10 the most fundamental of which is that such searches often find many more hypothetical crystal structures within the small energy range of possible polymorphism than known polymorphs. While relatively little progress has been made on modeling the kinetic factors that determine which thermodynamically feasible crystal structures will actually be observed, this field has inspired work on improving the thermodynamic modeling of organic crystal structures,11 their intermolecular forces,9,12 and the search methods.13 Because of their pharmaceutical and general importance, some confor-
10.1021/cg049922u CCC: $27.50 © 2004 American Chemical Society Published on Web 10/07/2004
1120
Crystal Growth & Design, Vol. 4, No. 6, 2004
mationally flexible molecules have been studied from the earliest days, but these have been acknowledged to be considerably more challenging. The simplest approach14 to conformationally flexible molecules is to use a force field that models both the inter- and the intramolecular forces. These have been extensively developed for modeling biological systems, and the crystal structure prediction work should benefit from the immense effort that has gone into force field development. Examples of crystal structure predictions on conformationally flexible molecules range from primidone15 to the hugely computationally demanding limited study on two diastereomeric salts of a chlorinesubstituted cyclic phosphoric acid and ephedrine.16 Initiating several searches with the molecule in different low energy conformers often ameliorates the problem of ensuring an adequate search. Hence, the main limitation on the accuracy of the results has often been identified as the accuracy of the force field. Crystal structure prediction can be particularly sensitive to whether the balance between the intermolecular forces and the intramolecular forces produces the correct molecular conformation and thus relative total energies. A recent study17 of the ability of variations of an atomatom force field to reproduce the crystal structures of a wide range of pharmaceutical type compounds found many cases where the unphysical distortion of the molecular conformation prevented the crystal structure being even qualitatively reproduced. An impressive series of crystal structure prediction studies18,19 on highly flexible alcohols and sugars also led to the conclusion that force field accuracy was the major issue. This led to the development20 of a highly accurate ab initio-based anisotropic, polarizable, atomatom model potential, derived from ab initio calculations on methanol and its dimers and trimers, for calculating the intermolecular lattice energy, Ulatt. This potential was seen to be transferable to other alcohols, alkanes, and ethers and gave superior results in crystal structure prediction studies21 to standard empirical force fields. Further efforts22 to improve the accuracy of the predictions for glycol and glycerol led to the use of this intermolecular potential with increasingly sophisticated intramolecular force fields. The estimates of the intramolecular energy penalty ∆Eintra incurred by the force field-estimated distortions of the molecule from its lowest energy conformation in response to the crystal packing forces22 were finally refined by ab initio calculations. However, in principle, the conformations should also be determined ab initio, as the most stable 0 K static crystal structure is given by optimizing ∆Etot ) Ulatt + ∆Eintra, by balancing the inter- and intramolecular forces. Further progress on glycol and glycerol23 has been made toward this end in using both ab initio conformational energies (∆Eintra) and associated forces, within the final optimization process. This use of highly accurate models for both the inter- and the intramolecular forces proved an enormous improvement on traditional force fields, but extension of this scheme requiring so many ab initio energies and derivatives to larger organic molecules is currently computationally prohibitive. Thus, in this paper, we assess the utility and limitations of an intermediate approach, based on high level
Ouvrard and Price
Figure 1. Atom definitions and torsion angles for aspirin. For the planar conformation illustrated, τ1 ) 180°, τ2 ) 0, and τ3 ) 180°.
ab initio determinations of the gas phase conformers and their relative energies as ∆Eintra. These are used in rigid body searches for the most stable crystal packings, using a lattice energy search with an accurate distributed multipole model to estimate the relative lattice energies. Clearly, the success of this approach depends on how close the molecular conformation in the solid state is to a gas phase conformation, and second, how sensitive the lattice energy search is to these packing-induced differences. The degree of conformational change induced by the packing forces is often modest24 although obviously very dependent on the strength of the torsional potential about the most flexible single bonds relative to the packing forces. Hence, as a preliminary and discussion aid, it is also necessary to consider the torsional potential energy surfaces. Aspirin has been chosen for this study as a simple molecule with three main torsion angles shown in Figure 1. Rotation of the hydroxyl and methyl groups was implicitly considered. There is significant dynamical motion of the methyl group in the crystal, with 〈φ2〉 ≈ 600°2 around room temperature, which can be wellmodeled without considering its intermolecular contacts.25 The carboxylic group is expected to be coplanar with the aromatic ring unless repelled by the orthosubstituent, whereas the ester group has more conformational flexibility. Aspirin is of particular interest, as a crystal structure prediction study26 using the DREIDING force field showed that the most stable computed crystal structures were for a planar conformation of aspirin but that the structure with the most stable lattice energy approximated the experimental structure. However, the authors noted that the force field used in the crystal structure prediction gave the planar conformer as the most stable in the gas phase, whereas semiempirical AM1 calculations gave the crystal structure conformation as more stable. Thus, the prediction26 that polymorphs of aspirin might be found if routes could be found to promote the stability of more planar conformers rests on the accuracy and balance of interand intramolecular energies. The relative stability of the gas phase conformers of aspirin has been extensively studied27 and discussed using Hartree-Fock (HF) and density functional ab initio methods, giving a variety of low energy conformers, none of which are planar. Moreover, the gas phase conformer that is most similar to that found in the crystal structure25,28,29 is about 3.5 kJ mol-1 less stable than the global minimum. Although there is significant elongation in the thermal ellipsoid of the hydroxyl
Conformationally Flexible Molecules of Aspirin
proton above 200 K,25 this is likely to reflect anharmonicity in the hydrogen bond potential rather than static or dynamic disorder in the carboxylic acid dimer. Hence, this recent neutron study implies that there is only the one conformation in the crystal structure.25 Thus, in this study, we seek to extend the gas phase conformational analysis and then explicitly consider the ability of rigid low energy conformers, both minima and planar, to form low energy crystal structures. Aspirin is a good candidate for crystal structure prediction studies as it has been heavily studied and frequently crystallized under a variety of conditions, because of its therapeutic and commercial importance. Although considerable variations in some physical properties, such as morphology and dissolution rates, are observed for aspirin, these do not constitute evidence for a new polymorph,30,31 and despite various efforts, only the one crystal structure has been established by diffraction. It therefore seems unlikely that there are any long-lived polymorphs. Materials and Methods To establish whether it was likely that Glaser27 had located all significant conformational energy minima, and also the nature of the potential wells around the observed conformer, an initial scan of the conformational potential energy surface for the acetyl substituent of aspirin (τ1 ) C1-C2-O3-C8 and τ2 ) C2-O3-C8-O4) was carried out at the AM1, PM3, and B3LYP/6-31G(d) levels of calculation, with all other degrees of freedom fixed. Because it is assumed that the carboxylic acid substituent does not deviate significantly out of the plane of the phenyl ring,32 the corresponding torsion angle τ3 ) C2C1-C7-O2 was held fixed at τ3 ) 0° and τ3 ) 180°. These analyses with both semiempirical and ab initio methods allowed comparison of the approximate positions of the minima in the conformational energy and analysis of the consistency between methods. Indeed, the previous work of Payne et al.26 showed that AM1 and a DREIDING force field approach led to very different low energy conformations. The minimum energy conformers of aspirin were optimized at the HF, density functional B3LYP, and Moller-Plesset MP2 levels of theory using the 6-31G(d,p) basis set, using Gaussian98.33 This extends the study of Glaser27 to contrast the treatment of electron correlation between the MP2 and a density functional method. Also, geometrical optimizations were carried out using a slightly larger basis set, with polarization functions on both heavy atoms and hydrogen, which should give a better description of the charge distribution in the hydrogen-bonding protons. Constrained optimizations of different planar conformations were also performed at the B3LYP level to test the suggestion26 that these transition state conformations were of sufficiently low energy that they might be produced by crystal packing forces. The crystal structure predictions used the B3LYP optimized molecular structures. In addition, for the comparison with the search with optimized structures, the molecular structure of the room temperature crystal28 was used directly, with only the hydrogen atom positions adjusted to elongate the bonds to the standard neutron-diffraction-based values34 of 1.08 Å for C-H and 1.02 Å for O-H. While genuine crystal structure predictions using rigid molecules have to use a computed molecular structure, a significant proportion of work on crystal structure prediction takes the rigid molecular structure from a room temperature crystal structure. To model the dominant electrostatic contribution to the lattice energy accurately, the molecular charge distribution of each conformer was represented by a set of atomic multipoles up to hexadecapole. These had been obtained by a distributed multipole analysis (DMA)35 of the B3LYP charge density, calculated using the program GDMA36 on the Gaussian9833 charge density file. The electrostatic contribution to
Crystal Growth & Design, Vol. 4, No. 6, 2004 1121 the lattice energy included all terms in the atom-atom multipole series up to R-5, with the charge-charge, chargedipole, and dipole-dipole terms being evaluated by Ewald summation. The remaining terms were calculated by direct summation up to a molecule-molecule separation of 15 Å. The empirical repulsion-dispersion potential has the form
U)
∑
Uik )
i∈1,k∈2
∑
(AιιAκκ)1/2 exp[-(Bιι + Bκκ)Rik/2] -
i∈1,k∈2
(CιιCκκ)1/2 Rik6
(1)
where atom i in molecule 1 is of type ι and atom k in molecule 2 is of type κ. The parameters for the C, HC, and O atoms were taken from the empirically fitted potential of Williams et al.,37,38 and the hydroxyl protons HO parameters were empirically fitted as an addition to this parameters set for use with a DMA electrostatic model.39 The crystal structure prediction searches used MOLPAK40 to generate hypothetical dense packings of the various rigid molecular structures. MOLPAK performs a systematic grid search on orientations of the central molecule in 29 common coordination geometries of organic molecules, belonging to the space groups P1, P1 h , P21, P21/c, P21212, P212121, Pna21, Pca21, Pbca, Pbcn, C2/c, Cc, and C2. Approximately the 50 densest packings in each coordination type were then used as starting points for lattice energy minimization, using the DMA-based model potential. All lattice energy minimizations were carried out using DMAREL3.02,41,42 the extension of DMAREL43 in which the minimizations are performed within the space group symmetry. After minimization within the symmetry constraints of the space group, the eigenvalues of the approximate second derivative matrix were examined to ensure that true minima had been found, and any transition states were discarded. The reduced cell parameters of the lattice energy minima were calculated using PLATON,44 and duplicate structures were discarded. Many of the searches with the different conformers were carried out as part of the development of a prototype Grid computational service.45 The service is implemented following the open grid services architecture. This service allows one submission to perform the complete MOLPAK search and DMAREL optimizations and gathering of the results for one conformation into a database as one procedure. The automatic distribution of the individual steps over the available processors on a Beowulf type cluster, when each minimization and MOLPAK packing type search takes a variable amount of CPU, produces an enormous savings in human time. This e-Science work will make searches with a variety of conformations and a computationally demanding model intermolecular potential feasible. These searches were carried out with those B3LYP ab initio conformations that were sufficiently low in energy that it was plausible that any loss in conformational energy ∆Eintra could be balanced by improvements in the lattice energy, Ulatt. All of the unique structures within 7 kJ mol-1 of the global minimum of ∆Etot ) Ulatt + ∆Eintra were reminimized from the initial MOLPAK starting point using DMAREL3.1, which uses accurate second derivatives and allows calculation of both the elastic constant matrix46 and the rigid body intermolecular harmonic phonon frequencies.47
Results The level of calculation of the potential energy surface scan, semiempirical or ab initio, does not seem to influence the approximate positions of the minima in a rigid model scan. Indeed, conformational maps of torsion angles τ1 and τ2 calculated at AM1, PM3, and B3LYP/ 6-31G(d) (an example corresponding to τ3 ) 180° is illustrated in Figure 2) levels have similar shapes and are also in good agreement with the potential energy surface established by Payne26 at the AM1 level of
1122
Crystal Growth & Design, Vol. 4, No. 6, 2004
Ouvrard and Price Table 1. Relative Energies and Dipole Moments of Aspirin Conformers as Determined by Ab Initio Optimization MP2/6-31G(d,p) B3LYP/6-31G(d,p) HF/6-31G(d,p) structurea 1a 2a planar A 4a planar B 5a planar C 1b 2b 3a 4b 3b
Figure 2. Map of variations in the B3LYP/6-31G(d) conformational energy of aspirin with angles τ1 and τ2 with all other bond lengths and angles held fixed as in the metastable conformer 2a, which is closest to the crystal structure conformer. The energy plotted is relative to the energy minimum of conformer 2a. Because τ3 ) 180°, the secondary minimum on this plot is conformer 2b. Relative conformational energies are color-coded in 5 kJ mol-1 intervals, from red squares denoting 0-5 kJ mol-1 to blue squares for 45-50 kJ mol-1. Conformational energies above 50 kJ mol-1 are gray-colored squares.
calculation. This shows that the set of conformational energy minima found by Glaser27 is almost certainly complete. The scanning map of torsion angles τ1 and τ2 calculated at the B3LYP/6-31G(d) level with the observed conformation of the carboxylic acid group (τ3 ) 180°, Figure 1) shows that a low energy conformation is found for τ1 ∼ (90° and τ2 ∼ 0°. This conformer 2a corresponds to the approximate orientation of the acetyl substituent observed in the experimental crystal structures. It also reveals that the acetyl is a quite flexible substituent as τ1 and τ2 can vary by about (20° about this minimum with less than 5 kJ mol-1 conformational energy penalty. Thus, the acetyl group could easily rotate quite significantly about both bonds C2-O3 and O3-C8 under the influence of crystal packing forces. A second minimum for the observed conformation of the carboxylic acid group (Figure 2) with s-cis conformation about bond O3-C8 (τ1 ∼ (90° and τ2 ∼ 180°, corresponding to conformer 2b) is found to be much higher in energy (Figure 2 and Table 1) than estimated by the DREIDING force field or AM1 calculations.26 All of the nine optimized gas phase conformers have the ester group significantly rotated about τ1 to be out of the plane of the aromatic ring. The minima that are s-cis about O3-C8 (conformations of type b in Table 1) are much higher in energy, with the proximity of the C8-C9H3 group forcing the carboxylic group out of the plane of the ring. The two lowest energy conformations differ mainly in the position of the s-trans carboxylic acid proton (Figure 3), but the possibility of an internal hydrogen bond for the s-cis conformation stabilizes some of these higher energy conformations. The B3LYP and HF energies and structures are in close agreement with those determined by Glaser,27 as expected from the
∆Eb (kJ mol-1)
µ (D)
0.00 2.86
2.18 3.44
15.70
5.50
17.17
5.55
20.96 22.52 23.90 31.06 60.23
6.37 5.04 4.53 7.07 9.02
∆Eb (kJ mol-1)
µ (D)
0.00 3.47 11.94 12.09 12.12 15.27 15.40 21.08 21.83 26.30 28.45 60.45
2.02 2.64 1.18 4.61 1.44 5.19 3.86 5.36 3.84 3.83 6.07 7.68
∆Eb (kJ mol-1) µ (D) 0.00 4.59
2.37 2.94
23.83
5.07
24.71
5.32
28.92 30.00 31.35 44.34 76.68
6.09 4.54 4.35 6.54 8.66
a The structures are denoted as in ref 27, with the a structures being s-trans about O3-C8 giving a methyl group sticking out, and the b structures s-cis with C8-O4 protruding. The three lowest energy minima and the planar transition states are shown in Figures 3 and 5, respectively. b The energy is given relative to the lowest energy gas phase structure 1a. Energies in italics are those of first-order transition states.
Figure 3. Molecular structures of the three lowest energy gas phase conformers optimized at the B3LYP/6-31G(d,p) level.
minor difference in the basis set. However, the MP2 results, while preserving the general ordering of the conformers, differ quite significantly from the density functional results both in energies (Table 1) and in conformation. Indeed, Figure 4 clearly shows that there are significant differences in the orientation of the ester and carboxylic acid groups relative to the aromatic ring between the MP2 and the B3LYP conformations of minimum 2a. The B3LYP structure is close to the experimental structures and, in particular, in very good agreement with the 20 K molecular structure, as might be expected as both approximate the 0 K static structures. The experimental room temperature X-ray structure28 has a C7-O1 bond length of 1.289 Å and a C7d O2 bond length of 1.239 Å, which differ somewhat from the 20 K neutron values25 of 1.309 and 1.223 Å, respectively. They also contrast with the B3LYP ab initio values of 1.356 and 1.215 Å or MP2 values of 1.361 and 1.220 Å for C7-O1 and C7dO2 bond lengths, respectively. Hence, the experimental room temperature molecular model differs from the ab initio optimized model in terms of both bond lengths and angles. The planar conformations of aspirin (Figure 5) are all found as first-order transition states. Their energies are about 12 kJ mol-1 above the global minimum (Table 1), confirming the suggestion made by the semiempirical calculations,26 that the DREIDING force field was incorrect in predicting a stable planar conformer. However, because planar molecules may have an advantage in producing dense, energetically favorable, crystal
Conformationally Flexible Molecules of Aspirin
Crystal Growth & Design, Vol. 4, No. 6, 2004 1123
Figure 5. Structures of the three lowest energy gas phase planar structures optimized at the B3LYP/6-31G(d,p) level. These structures are all first-order transition states.
Figure 6. Reproduction of the experimental crystal structure of aspirin by the computational model. Two views overlaying the room temperature crystal structure28 in red, with the lattice energy minimum obtained with the DMA-based intermolecular potential using the B3LYP/6-31G(d,p) molecular conformation 2a (black). The black figure is also the closest approximation to the experimental crystal structure found in the search with the ab initio optimized conformer 2a.
Figure 4. Superposition of the crystal structure conformations at 20 K (green) and room temperature (black) with the gas phase conformer 2a optimized at B3LYP/6-31G(d,p) (red) and MP2/6-31G(d,p) (blue) levels of theory. The dihedral angles between the mean plane of the aromatic ring and the carboxylic acid (C7-O1-O2) are 1.7, 2.5, 5.1, and 21.6°, and the dihedral angles between the mean plane of the aromatic ring and the acetyl group (C2-O3-C8-C9) are 81.2, 84.5, 77.9, and 63.8°, respectively.
packings, the lowest energy planar conformations A and B were used in a crystal structure prediction calculation. Crystal Structure Reproduction. The reproduction of the known crystal structure by lattice energy minimization (Table 2) confirms the adequacy of the DMA plus empirical exp-6 potential for this study. This potential, a variation with different HO parameters,48 plus the DREIDING force field used in the previous crystal structure prediction study26 and a simple rigid body exp-6 plus point charge model28 all give reasonable
reproductions of the crystal structure considering the inherent errors in comparing lattice energy minima with room temperature crystal structures.49 Nevertheless, the use of the nearest ab initio optimized molecular structure (conformer 2a) does lead to a slightly poorer reproduction of the room temperature crystal structure (ExptMinOpt2a, Table 2 and Figure 6) than when the molecular structure taken from the crystal structure28 is used (ExptMinExpt, Table 2). The use of the gas phase molecular structure gives a less favorable lattice energy, which is consistent with the small conformational differences being induced by the crystal packing forces. It is notable that despite the similar optimized crystal structures, the three different intermolecular potential models give very different lattice energies. The DREIDING force field gives by far the largest lattice energy but has previously been observed to tend to overestimate lattice energies for many hydrogen-bonded crystal structures.50 Because the electrostatic contribution is about half of the total lattice energy, the most accurate representation of the charge distribution used in the DMA model suggests that this value should be more accurate.
Table 2. Reproduction of the P21/c Crystal Structure of Aspirin by Lattice Energy Minimization experimental 20 K neutron25 a (Å) b (Å) c (Å) β (°) RMS (%)a V (Å3) Ulatt (kJ mol-1) a
11.186 6.540 11.217 96.07 816.0
literature computational
room temp X-ray28
6-exp + H-bond28
11.430 6.591 11.395 95.68 1.4 854.2
11.434 6.390 11.169 97.70 1.8 808.7 -122.8
prediction DREIDING26 11.42 6.69 11.22 95.04 1.6 -184.9
RMS: percentage error in the lattice constants, relative to the 20 K structure.
current potential + rigid conformer ExptMinExpt
ExptMinOpt2a
11.541 6.582 11.307 95.31 1.7 855.3 -114.9
11.389 6.758 11.350 95.86 2.0 869.0 -106.1
1124
Crystal Growth & Design, Vol. 4, No. 6, 2004
Figure 7. Hypothetical crystal structures generated using different gas phase conformations in the initial MOLPAK/ DMAREL search. The energy Etot ) Ulatt + ∆Eintra is plotted against the cell volume per molecule, where the intramolecular energy contribution ∆Eintra is the B3LYP/6-31G(d,p) conformational energy difference (Table 2) for each of the ab initio optimized minima 2a (orange circles) and 1a (maroon circles) and planar transition states planar A (blue circles) and planar B (green circles).
Figure 8. Lowest energy hypothetical crystal structures generated using different ab initio gas phase conformations in the initial MOLPAK/DMAREL search (i.e., detail from Figure 7), denoted by molecular conformer and space group [2a, P21/c ([); 2a, C2/c (9); 1a, P21/c (2); and 1a, P1 h (b)].
Crystal Structure Prediction Search. The total packing energies (∆Etot ) Ulatt + ∆Eintra) of the hypothetical structures generated in the search with the low energy optimized conformations are depicted in Figure 7 (detail of the area around the global minimum is given in Figure 8). It shows that the most stable conformer in the gas phase, 1a, cannot achieve such a low lattice energy as conformer 2a does, so that both conformers produce crystal structures with energies around the global minimum in ∆Etot. The planar structures can generate crystal structures, which have very similar
Ouvrard and Price
Figure 9. Calculated X-ray powder diffraction patterns for the room temperature experimental structure28 (green), the closest structure found in the search using the ab initio optimized conformer 2a (AI11 and ExptMinOpt) (orange) and another hypothetical structure, which is indistinguishable on energy but has some minor differences in the packing (2a_AK22) (blue).
lowest lattice energies to those generated by conformer 2a, but their intramolecular energy penalties make a crystal structure based on a planar conformation energetically unfavorable. The next low energy gas phase conformer, 4a, is considerably less stable and yet is stabilized by an intramolecular hydrogen bond resulting in less opportunity for intermolecular hydrogen bonding to occur in the crystal. Hence, it and the other much higher energy conformations were not used in the crystal structure prediction search. Figure 8 shows that there are three unique crystal structures close to the global minimum in ∆Etot, separated by 2 kJ mol-1 from the majority of other crystal structures. One of these low energy structures (2a_AI11, Table 3) is the experimental crystal structure, as reproduced as accurately as the computational model allows (i.e., ExptMinOpt2a in Table 2). At about the same total packing energy (assuming the B3LYP value for ∆Eintra), the most stable conformer of aspirin, 1a, gives rise to a crystal structure (1a_AM6, Table 3) in which there are hydrogen bonds from the carboxyl proton to the acetyl-carbonyl group. Both the experimental and the lowest energy crystal structure (2a_AK22, Table 3) have the carboxylic acid dimer as the hydrogenbonding motif. However, although the global minimum structure (2a_AK22) has a similar molecular arrangement to the marginally less stable experimental crystal structure, it has a low shear elastic constant, Cij, of 0.15 GPa, implying that it is so readily deformed that there may be problems in its growth. The similarities between the two structures, also apparent in the calculated X-ray powder diffraction patterns (Figure 9), suggest that the hypothetical structure (2a_AK22) might readily transform to the experimental structure if it were thermodynamically more stable. Hence, it is worth investigating the sensitivity of the relative energies to the molecular model. Comparison of Crystal Structure Prediction Using Ab Initio Optimized Conformers and Experimental Molecular Structures. The low energy crystal structures generated in a search with the experimental molecular structure are tabulated in Table
Conformationally Flexible Molecules of Aspirin
Crystal Growth & Design, Vol. 4, No. 6, 2004 1125
Table 3. Low Energy Crystal Structures Found in the Search with Ab Initio Optimized Molecular Structuresa cell parameters structureb
space Ulatt (kJ Etotc (kJ density group mol-1) mol-1) (g cm-3)
a (Å)
b (Å)
c (Å)
hydrogen bonds d(O‚‚‚O) (Å) and (O‚‚‚HO) (°) lowestd R (°) β (°) γ (°) Cij (GPa) O4‚‚‚H8-O1 O2‚‚‚H8-O1 motif
2a AK22 P21/c -106.3 1a AM6 P21/c -102.7
-102.9 -102.7
1.381 1.428
12.124 6.696 11.475 90.0 68.4 90.0 9.800 9.152 9.434 90.0 82.0 90.0
0.15 4.72
2.868 (153)
P21/c P21/c P21/c C2/c
-106.1 -100.8 -100.7 -104.0
-102.7 -100.8 -100.7 -100.5
1.377 1.405 1.329 1.418
11.388 9.828 7.169 17.649
90.0 90.0 90.0 90.0
4.52 5.40 2.25 2.55
2.877 (132)
1a AB25 P1 h 1a AK7 P21/c
-100.3 -100.1
-100.3 -100.1
1.349 1.350
6.094 6.982 10.850 79.6 78.6 83.3 11.052 6.829 11.835 90.0 84.1 90.0
1.86 3.02
2a 1a 1a 2a
AI11 FC2 AK15 DC46
6.758 6.993 7.172 8.953
11.350 13.375 17.538 15.199
90.0 90.0 90.0 90.0
95.9 67.9 93.2 44.7
2.852 (170)
dimer infinite chains 2.857 (169) dimer 2.879 (155) dimer 2.837 (170) dimer 2.889 (118) three centers 2.817 (172) dimer 2.848 (167) dimer
a The structure in bold (2a_AI11) corresponds to the experimental crystal structure (ExpMinOpt, Table 2), and the structure in italics seems unlikely to be an observable polymorph as it has a very small shear elastic constant. b The crystal structure is denoted by the molecular conformer and one of the MOLPAK-generated starting structures that resulted in this minimum. c The sum of the lattice energy and the B3LYP/6-31G(d,p) relative conformational energy from Table 1. d The smallest eigenvalue of the shear submatrix of the elastic constant tensor.
Table 4. Low Energy Crystal Structures Found in the Search with the Experimental Molecular Structurea experimental molecular structure space group P21/c P21/c P21/c P1 h P21/c P21/c P1 h
ab initio conformation 2a
structure
Ulatt (kJ mol-1)
cell vol/ Z (Å3)
a (Å)
b (Å)
c (Å)
R (°)
β (°)
γ (°)
structure
Ulatt (kJ mol-1)
AI11 AI12 FC49 AB33 AK44 AK24 CA8
-114.9 -114.9 -108.8 -108.7 -108.6 -108.3 -108.1
213.8 214.0 215.1 220.2 225.1 226.6 221.0
11.54 12.32 9.01 6.59 7.31 11.66 8.83
6.58 6.52 8.83 5.99 6.91 6.84 6.75
11.31 11.47 10.83 12.00 17.85 17.88 7.44
90.00 90.00 90.00 110.00 90.00 90.00 87.04
95.31 68.30 93.00 86.11 90.53 39.50 90.22
90.00 90.00 90.00 98.48 90.00 90.00 87.13
AI11 AK22 FC43 AB43 AK15 AK19 CA34
-106.1 -106.3 -101.4 -101.4 -101.9 -100.7 -100.6
a The structure in bold (expt_AI11) corresponds to the experimental structure (ExptMinExpt, Table 2). The final columns are the corresponding structures found in the search with the ab initio conformer 2a (mainly in Table 3), which has an identical hydrogen bond motif and a similar arrangement of the molecules with very similar cell parameters.
4. It shows that using the experimental molecular structure, the experimental crystal structure is found as the most stable (i.e., expt_AI11 matches ExptMinExpt). It is notable that there is a structural correspondence between most of the low energy structures found with the experimental molecular conformation and those found using the ab initio conformer 2a (Table 4). However, the small differences (Figure 4) in conformation between the experimental molecular structure and the ab initio optimized conformation 2a produce a significant destabilization of the lattice energies, of the order of 8 kJ mol-1, as already seen for the experimental crystal structure in Table 2. The difference between the two lowest energy structures and the other less stable structures increases from 5 to 6 kJ mol-1; that is, the experimental structure increases this differentiation, and there are small relative energy shifts within these general trends. In particular, the experimental crystal structure is now found as the most stable structure in lattice energy. This is in agreement with the PROMET search, performed with the experimental molecular structure and a very different intermolecular potential, which found the experimental structure as the global minimum.51 It is worth noting that all of the structures in Table 4 are based on the carboxylic acid dimer motif, which may account for the relative energy shifts being so small. Discussion The MOLPAK/DMAREL searches with ab initio optimized molecular structures give a clear prediction of the compensation between the contributions to the total
crystal packing energy (∆Etot ) Ulatt + ∆Eintra). A less stable molecular conformation can be compensated for by increased lattice energy and therefore can be observed in a crystal structure. This results in the experimental structure being found as one of the three most stable structures according to the total packing energy (∆Etot ) Ulatt + ∆Eintra). However, even the small differences between the gas phase local minimum structure 2a and X-ray molecular crystal structure are sufficient to affect the relative stability of different crystal structures. The relative energies of the most stable crystal structures for the lowest energy and observed conformers of aspirin will be affected by errors in the intermolecular potential because they adopt very different hydrogen-bonding motifs. The energy ordering will be further affected by the errors in the ab initio relative energies, ∆Eintra, let alone the neglect of temperature effects. Because the success of this study is so dependent on small energy differences, little confidence can be placed in the relative orderings of Etot for the three structures around the global minimum in Etot as a guide to thermodynamic stability. There may be kinetic reasons favoring the formation of the carboxyl dimers in the crystal structure, over the alternative hydrogen-bonded chain structures, since the phenomenon of polymorphism implies that thermodynamic stability is not the sole factor determining which crystal structures are observed. The rotation of the carboxyl group (τ3), to give the less stable conformer 2a, allows the formation of carboxylic acid dimer crystal structures, which are significantly more stable. Indeed, the lowest energy structure based on the carboxylic dimer of the 1a conformer, 1a_FC2 in Table 3, estimates
1126
Crystal Growth & Design, Vol. 4, No. 6, 2004
this energy difference as at least 5 kJ mol-1. If the growth units of the crystal are carboxylic acid dimers, then the occurrence of the less stable conformer is likely. Overall, if one had to select a few crystal structures from those produced by the search with the ab initio conformations (i.e., those in Table 3) and make a prediction, as in the blind tests,5,6 the crystal structure of aspirin might well have been successfully predicted. However, this prediction could not be made with much confidence. Although we do rule out the possibility of a polymorph with a planar conformation of aspirin, other polymorphs based on variations of conformers 1a or 2a and using either the carboxylic acid dimer or the acetyl group hydrogen-bonding motifs cannot be confidently excluded. Hence, this approach to crystal structure prediction based on lattice energy minimization using all of the low energy gas phase conformations is worth trying for molecules that can adopt a few distinct low energy conformations. The key criterion for suitability is whether the packing forces could not substantially modify the ab initio conformations. A comparison of 12 molecular substructures in the crystalline state against HF/STO3G conformational potential energy curves for model compounds24 showed that torsion angles with a significant strain energy (>1 kcal mol-1) appear to be very unusual in crystal structures. The average observed torsion angle was generally within a few degrees of the ab initio optimized angle for the functional groups considered, although individual molecules often had much larger deviations. In a test of whether the discrepancies between the sublimation energy and the lattice energies calculated by some force fields could be attributed to a conformational energy change, the differences between the crystal and the HF/6-31G(d,p) conformations were found to be small.50 More generally, comments have often been made on the similarity between the ab initio prediction and the observed conformation, as for example the three conformational polymorphs of di-µ-chlorotetrakis(1-methylboratabenze)diyttrium.52 Hence, in such cases, it would be worth trying a crystal structure prediction study based on the gas phase low energy conformers, although the success will be critically dependent on the sensitivity of the relative lattice energies to the errors in the conformation. The chances for a successful prediction are best when the differences in the energies of the low energy crystal structures are large as compared with the changes in those energies arising from variations in the molecular structure. Nevertheless, this approach cannot work for cases where the crystal packing obviously determines the molecular conformation, for example, when an internal hydrogen bond occurs in the gas phase structure, and the proton donor group can distort to form an intermolecular hydrogen bond in the crystal (e.g., o-acetamidobenzamide53). Discussions of the differences between crystal and gas phase conformations generally assume that the ab initio method gives a very good approximation to the gas phase structure. It is a matter for concern that the conformational differences between the MP2 and the B3LYP optimized structures are far more significant than the differences between the B3LYP and the experimental structures. Both methods give a far better prediction of the relative position of the acetyl group than HF theory, which can be attributed to both
Ouvrard and Price
correlated methods representing the nonbonded dispersion interaction. The agreement of the B3LYP and low temperature crystal structure conformations of aspirin could have resulted from some fortuitous cancellation of errors. It is not clear that the MP2 structure should be a worse approximation to the actual gas phase structure, although it may be more affected by basis set superposition error.54 Thus, there seems to be more uncertainty in the gas phase structure from the ab initio theory than packing induced differences in conformation. However, given the shallowness of the potential energy surface in this region (Figure 2) and the different dynamical motions of the molecules in both the gas and the crystalline phases, the differences in conformation are perhaps only a practical problem for this approach to crystal structure prediction. The considerable sensitivity of the lattice energy to the conformation shows that a large number of low energy conformations should be used in the search. Indeed, the conformational energy scan (e.g., Figure 2) suggests that a search with many other low energy conformations, for example, variations on conformer 2a with τ1 and τ2 changed by say 10-20°, would need to be performed for a genuine prediction. Because the experimental molecular conformation of aspirin was closer to the ab initio minimum (Figure 4) than such grid points, these variations have not been performed. A similar grid search would have been necessary to consider rotations of the carboxylic acid group had the internally hydrogen-bonded conformer 4a been much closer to a gas phase minimum. Thus, ab initio conformers can only be used as rigid probes in crystal structure prediction, with considerable care, and consideration of whether there are readily distortable hydrogen-bonding groups. Conclusions The combination of rigid body molecular models and lattice energy minimization is certainly successful for determining which conformations are more likely to have low energy crystal structures. In the case of aspirin, a fairly successful crystal structure prediction was possible, because the energy gap between the few lowest energy structures and the other hypothetical structures was larger than the changes in relative lattice energies produced by replacing the ab initio optimized conformation with the observed molecular structure. Unfortunately, it appears that quite high levels of ab initio optimization are needed for determining the range of low energy conformations. This, combined with the need for accurate intermolecular potentials for the lattice energy, implies considerable computational costs, and inevitably, there will be compromises between these factors and the completeness of the search in terms of the range of conformations and crystal packing types. Ideally, crystal structure prediction of conformationally flexible molecules requires a combined model where the inter- and intramolecular energies are calculated very accurately and the forces are balanced in the optimization. This study generally confirms the conclusion derived from the glycol and glycerol studies11 that high level ab initio conformational forces and energies will need to be combined with the best quality intermolecular force fields to obtain sufficiently accurate Etot energies in many cases. However, the required accuracy
Conformationally Flexible Molecules of Aspirin
is not an absolute but is relative to the energy differences between the hypothetical structures. This study does suggest that for larger molecules, more typical of pharmaceuticals, there may only be a few conformations that are both low in energy and can pack effectively. Also, because energy alone is not a predictor of crystal structures, a method that effectively contrasts the ability of different conformations to pack effectively into crystal structures, even if there is considerable uncertainty in the relative Etot, is worth pursing toward an understanding of conformational polymorphism. Acknowledgment. Particular thanks to Ben Butchart and Wolfgang Emmerich of UCL Computer Science Department for the development of the OGSA service to run MOLPAK/DMAREL crystal structure prediction searches. This work was funded by the EPSRC project “e-Science Technologies in the Simulation of Complex Materials”. References (1) Bernstein, J. Polymorphism in Molecular Crystals; Clarendon Press: Oxford, 2002. (2) Yu, L.; Stephenson, G. A.; Mitchell, C. A.; Bunnell, C. A.; Snorek, S. V.; Bowyer, J. J.; Borchardt, T. B.; Stowell, J. G.; Byrn, S. R. J. Am. Chem. Soc. 2000, 122, 585-591. (3) Chemburkar, S. R.; Bauer, J.; Deming, K.; Spiwek, H.; Patel, K.; Morris, J.; Henry, R.; Spanton, S.; Dziki, W.; Porter, W.; Quick, J.; Bauer, P.; Donaubauer, J.; Narayanan, B. A.; Soldani, M.; Riley, D.; McFarland, K. Org. Process Res. Dev. 2000, 4, 413-417. (4) Bauer, J.; Spanton, S.; Henry, R.; Quick, J.; Dziki, W.; Porter, W.; Morris, J. Pharm. Res. 2001, 18, 859-866. (5) Lommerse, J. P. M.; Motherwell, W. D. S.; Ammon, H. L.; Dunitz, J. D.; Gavezzotti, A.; Hofmann, D. W. M.; Leusen, F. J. J.; Mooij, W. T. M.; Price, S. L.; Schweizer, B.; Schmidt, M. U.; van Eijck, B. P.; Verwer, P.; Williams, D. E. Acta Crystallogr. B 2000, 56, 697-714. (6) Motherwell, W. D. S.; Ammon, H. L.; Dunitz, J. D.; Dzyabchenko, A.; Erk, P.; Gavezzotti, A.; Hofmann, D. W. M.; Leusen, F. J. J.; Lommerse, J. P. M.; Mooij, W. T. M.; Price, S. L.; Scheraga, H.; Schweizer, B.; Schmidt, M. U.; van Eijck, B. P.; Verwer, P.; Williams, D. E. Acta Crystallogr. B 2002, B58, 647-661. (7) Dunitz, J. D. Chem. Commun. 2003, 545-548. (8) Beyer, T.; Lewis, T.; Price, S. L. CrystEngComm 2001, 3, 178-212. (9) Gavezzotti, A. CrystEngComm 2002, 4, 343-347. (10) Gavezzotti, A. Modell. Simul. Mater. Sci. Eng. 2002, 10, R1R29. (11) van Eijck, B. P. J. Comput. Chem. 2001, 22, 816-826. (12) Gavezzotti, A. J. Phys. Chem. B 2003, 107, 2344-2353. (13) Bazterra, V. E.; Ferraro, M. B.; Facelli, J. C. J. Chem. Phys. 2002, 116, 5984-5991. (14) Verwer, P.; Leusen, F. J. J. Computer Simulation to Predict Possible Crystal Polymorphs; Verwer, P., Leusen, F. J. J., Eds.; Wiley-VCH: New York, 1998; Vol. 12, pp 327-365. (15) Payne, R. S.; Roberts, R. J.; Rowe, R. C.; Docherty, R. Int. J. Pharm. 1999, 177, 231-245. (16) Leusen, F. J. J. Cryst. Growth Des. 2003, 3, 189-192. (17) Brodersen, S.; Wilke, S.; Leusen, F. J. J.; Engel, G. Phys. Chem. Chem. Phys. 2003, 5, 4923-4931. (18) van Eijck, B. P.; Mooij, W. T. M.; Kroon, J. Acta Crystallogr. B 1995, 51, 99-103. (19) van Eijck, B. P.; Kroon, J. J. Comput. Chem. 1999, 20, 799812. (20) Mooij, W. T. M.; van Duijneveldt, F. B.; van Duijneveldtvan de Rijdt, J.; van Eijck, B. P. J. Phys. Chem. A 1999, 103, 9872-9882. (21) Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Phys. Chem. A 1999, 103, 9883-9890. (22) Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Am. Chem. Soc. 2000, 122, 3500-3505. (23) van Eijck, B. P. J. Comput. Chem. 2001, 22, 805-815.
Crystal Growth & Design, Vol. 4, No. 6, 2004 1127 (24) Allen, F. H.; Harris, S. E.; Taylor, R. J. Comput.-Aided Mol. Des. 1996, 10, 247-254. (25) Wilson, C. C. New J. Chem. 2002, 26, 1733-1739. (26) Payne, R. S.; Rowe, R. C.; Roberts, R. J.; Charlton, M. H.; Docherty, R. J. Comput. Chem. 1999, 20, 262-273. (27) Glaser, R. J. Org. Chem. 2001, 66, 771-779. (28) Kim, Y.; Machida, K.; Taga, T.; Osaki, K. Chem. Pharm. Bull. 1985, 33, 2641-2647. (29) Wilson, C. C. Chem. Phys. Lett. 2001, 335, 57-63. (30) Pfeiffer, R. R. J. Pharm. Pharmacol. 1970, 23, 75-76. (31) Mitchell, A. G.; Milaire, B. L.; Saville, D. J.; Griffiths, R. V. J. Pharm. Pharmacol. 1971, 23, 534-535. (32) Leiserowitz, L. Acta Crystallogr. B 1976, 32, 775-802. (33) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Chesseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayalla, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. GAUSSIAN 98, A6 ed.; Gaussian Inc.: Pittsburgh, 1998. (34) Allen, F. H.; Kennard, O.; Watson, D. G.; Brammer, L.; Orpen, A. G.; Taylor, R. J. Chem. Soc., Perkin Trans. 2 1987, S1-S19. (35) Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047-1064. (36) Stone, A. J. GDMA: A Program for Performing Distributed Multipole Analysis of Wave Functions Calculated Using the Gaussian Program System, 1.0 ed.; Stone, A. J., Ed.; University of Cambridge: United Kingdom, 1999. (37) Williams, D. E.; Cox, S. R. Acta Crystallogr. B 1984, 40, 404-417. (38) Cox, S. R.; Hsu, L. Y.; Williams, D. E. Acta Crystallogr. A 1981, 37, 293-301. (39) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J. Phys. Chem. 1996, 100, 7352-7360. (40) Holden, J. R.; Du, Z. Y.; Ammon, H. L. J. Comput. Chem. 1993, 14, 422-437. (41) Price, S. L.; Willock, D. J.; Leslie, M.; Day, G. M. DMAREL 3.02; http://www.ucl.ac.uk/∼ucca17p/dmarelmanual/ dmarel.html, Collaborative Computational Project Number 5, Daresbury Laboratory, Daresbury, United Kingdom, 2001. (42) Mitchell, J. B. O.; Price, S. L.; Leslie, M.; Buttar, D.; Roberts, R. J. J. Phys. Chem. A 2001, 105, 9961-9971. (43) Willock, D. J.; Price, S. L.; Leslie, M.; Catlow, C. R. A. J. Comput. Chem. 1995, 16, 628-647. (44) Spek, A. L. PLATON, A Multipurpose Crystallographic Tool; University of Utrecht: Utrecht, 2002. (45) Butchart, B.; Chapman, C.; Emmerich, W. Proceedings of the UK e-Science All Hands Meeting 2003; Cox, S. J., Ed.; EPSRC: United Kingdom, 2003. (46) Day, G. M.; Price, S. L.; Leslie, M. Cryst. Growth Des. 2001, 1, 13-26. (47) Day, G. M.; Price, S. L.; Leslie, M. J. Phys. Chem. B 2003, 107, 10919-10933. (48) Beyer, T.; Price, S. L. J. Phys. Chem. B 2000, 104, 26472655. (49) Beyer, T.; Price, S. L. CrystEngComm 2000, 2, 183-190. (50) Osborn, J. C.; York, P. J. Mol. Struct. 1999, 474, 43-47. (51) Gavezzotti, A.; Filippini, G. J. Am. Chem. Soc. 1995, 117, 12299-12305. (52) Zheng, X. L.; Wang, B.; Englert, U.; Herberich, G. E. Inorg. Chem. 2001, 40, 3117-3123. (53) Buttar, D.; Charlton, M. H.; Docherty, R.; Starbuck, J. J. Chem. Soc., Perkin Trans. 2 1998, 763-772. (54) Tuma, C.; Boese, A. D.; Handy, N. C. Phys. Chem. Chem. Phys. 1999, 1, 3939-3947.
CG049922U
