1. Introduction

Crystal structure prediction (CSP) is of considerable interest in the development and manufacture of organic solid-state materials (Bernstein, 2002), not least in the pharmaceutical industry (Price, 2004a), where the prediction of a more thermodynamically stable polymorph of a drug candidate during the early stages of development may avert expensive problems later in the process. The discovery of a more stable polymorph during production may result in the need to reformulate the drug before it is possible to continue marketing it, as happened for the anti-HIV drug Norvir (ritonavir; Chemburkar et al., 2000) from Abbott Laboratories, where the new form was a conformational polymorph in which the ritonavir molecule adopts an alternative geometry (Bauer et al., 2001).

There has been considerable progress in the ability to predict the crystal structures of rigid molecules by searching for the global minimum in the lattice energy (Day, Motherwell & Jones, 2005) even under blind test conditions (Lommerse et al., 2000; Motherwell et al., 2002; Day, Motherwell, Ammon et al., 2005). However, this process is often based on calculating the molecular structure by ab initio optimization (i.e. estimating the gas-phase conformation) and assuming that this conformation is unchanged in the crystal structure. Molecules that have the possibility of adopting a range of conformations may, however, adopt a suboptimal gas-phase conformer such that the intermolecular interactions in the crystal are improved. A rather extreme example, which is investigated in this study, is when an intramolecular hydrogen bond is present in the ab initio gas-phase global minimum energy conformer, but not in stable solid-state packing arrangements where the hydrogen-bond donor/acceptor pair are only involved in lattice-stabilizing intermolecular interactions.

Thus, the extension of CSP methods to conformationally flexible molecules, which is important because industrially useful molecules often have some degree of flexibility, involves two major additional challenges. Firstly, in addition to modelling the intermolecular interactions in the crystal sufficiently to give the intermolecular lattice energy, U_latt, to a high degree of accuracy, which remains very challenging (Gavezzotti, 2002; Price, 2004b), it is also necessary to evaluate the relative energy penalty for the intramolecular distortion of the molecule from the gas-phase optimized conformation, E_intra, and use (1) to compare the relative thermodynamic stability of the hypothetical crystal structures

Secondly, the dimensionality of the search is increased, with the need to consider explicitly the intramolecular degrees of freedom, as well as the space group, cell parameters, and molecular position and orientation. Thus, it is not surprising that, for the three flexible molecules included in the CSP blind tests so far there has been only one successful crystal structure prediction by one of the many groups taking part (Lommerse et al., 2000).

The crystal structure prediction of a flexible molecule can be approached using an empirical atom-atom force field to model both intra- and intermolecular forces, and optimizing the energy E_tot with respect to the atom coordinates within the crystal. This approach has been successful in some cases (Verwer & Leusen, 1998), although the barriers to significant conformational change within a crystal and the multitude of local minima make it very challenging to ensure that sufficient conformations are sampled (Leusen, 2003). However, there may not be a sufficiently accurate atom-atom force field available for the molecule of interest; for example, a study of the crystal structures of many pharmaceutically relevant molecules found that the force field distorted the conformation of some of the flexible molecules sufficiently that the energy minimum was qualitatively different from the experimental structure (Brodersen et al., 2003). A series of crystal structure prediction studies on sugars and alcohols (van Eijck et al., 1995; van Eijck & Kroon, 1999) found that conventional force fields were inadequate. These authors progressed (van Eijck et al., 2001), in the cases of glycol and glycerol, by simultaneously optimizing U_latt and E_intra in the final minimization, where U_latt was calculated from a high-quality ab initio-based model intermolecular potential, and E_intra was calculated ab initio for the isolated molecule in the conformation induced by the intermolecular forces. The huge computer resources required to perform adequate-quality ab initio calculations for larger molecules within a crystal structure optimization imply that this procedure could, at best, only be carried out for the final refinement of low-energy structures in the foreseeable future.

Nevertheless, it is clear that good-quality ab initio estimates of the intramolecular contribution to the energy will often be required. A force-field-based study of aspirin found that the most stable calculated crystal structures were for the computed planar gas-phase conformer (Payne et al., 1999). A later study with various ab initio methods showed that the planar conformations were transition states between two low-energy non-planar conformations. A crystal structure prediction based on the two DFT (density functional theory) optimized low-energy conformers and planar transition states did rationalize why crystalline aspirin adopts the non-planar conformation, which is only a secondary minimum in the gas phase (Ouvrard & Price, 2004).

In this work, a strategy for the crystal structure prediction of flexible molecules is explored, based on the use of a strategic set of high-quality ab initio rigid conformers (with measures of relative E_intra) to search for low-energy crystal structures (with associated U_latt). The key torsion angles of the molecule are identified, and the molecular potential energy surface is scanned as a function of these angles to identify the ranges of conformers that might plausibly occur in the crystal, because a particularly good intermolecular lattice energy could compensate for a sufficiently small decrease in molecular stability (increase in E_intra). This region of conformational space is coarsely represented using sets of values for the key torsion angles. Each of these conformers is used to construct about 400 close-packed crystal structures in the most common coordination types (an exploratory search), which are then optimized by minimizing U_latt using an accurate intermolecular potential. The sum of the separately evaluated ab initio conformational energy E_intra and the minimized lattice energy U_latt, (1), is used to identify which conformers can give rise to crystal structures with favourable values of E_tot. These structures will generally correspond to structures in which the conformers can form intermolecular hydrogen bonds in a close-packed structure, though the coarse definition of the torsion angles and the use of a rigid conformer make it likely that a more stable structure would result if the molecular conformer was further optimized. Thus, the low-energy crystal structures from the initial searches are then used to define promising conformers for subsequent exploratory searches using crystallographic insight to choose related conformers that are expected to improve the crystal packing, for example, to bias the geometry of the conformer so a particular intermolecular hydrogen bond can shorten. Finally, full searches over more space groups are performed on the conformers that give the more stable structures in the exploratory searches.

Piracetam (2-oxo-1-pyrrolidineacetamide, C₆N₂O₂H₁₀, see Fig. 1) is an ideal choice for such a study because it has conformational polymorphs in which the molecular structures are very different to the gas-phase molecular structure with its internal hydrogen bond. Piracetam and piracetam-like compounds are of interest to the pharmaceutical industry because of their cognition-enhancing ability (Altomare et al., 1995; Gualtieri et al., 2002), for example in relation to the treatment of Alzheimer's disease (Evans et al., 2004) and acute stroke (Ricci et al., 2000), and have been studied extensively, resulting in a wealth of experimental evidence for assessing the calculations. Three polymorphs of piracetam had crystal structures in the Cambridge Structural Database (Allen, 2002) at the start of this study (see Table 1). Despite some indication of there being a further three phases with high degrees of metastability obtained from the melt (Kuhnert-Brandstätter et al., 1994), these were not found in a manual solvent evaporation screen (Keats, 2001). Form I is metastable; it exists at high temperature and spontaneously transforms to form II within a few hours at room temperature. The structure was first solved from powder diffraction data (Louër et al., 1995) and later a 150 K single-crystal structure determination showed significant disorder of the C4 atom in the pyrrolidine ring (Fabbiani et al., 2005). Forms II and III are stable at lower temperature, with form II being most stable at room temperature (Céolin et al., 1996), although recent work shows that forms II and III are very similar in stability around room temperature (Blagden & De Matos, 2005).

Table 1 Experimental X-ray determinations of piracetam All determinations were performed at room temperature using single-crystal data, except BISMEV03, which was determined using powder data, and FabbianiI, which was determined at 150 K. BISMEV02 has a better R factor than BISMEV01 (both are form III). FabbianiI and FabbianiIV are new determinations that are not yet in the CSD. Form I I II III III IV Refcode BISMEV03 FabbianiI BISMEV BISMEV01 BISMEV02 FabbianiIV Reference Louër et al. (1995) Fabbiani et al. (2005) Admiraal et al. (1982) Admiraal et al. (1982) Galdecki & Glowka (1983) Fabbiani et al. (2005) Space group P21/n P21/n P21/n P21/n P21/c a (Å) 6.747 (2) 6.7254 (2) 6.403 (3) 6.525 (2) 16.403 (3) 8.9537 (11) b (Å) 13.418 (3) 13.2572 (4) 6.618 (4) 6.440 (2) 6.417 (1) 5.4541 (6) c (Å) 8.090 (2) 8.0529 (2) 8.556 (6) 16.463 (5) 6.504 (1) 13.610 (4) (°) 90 90 79.85 (3) 90 90 90 (°) 99.01 (3) 98.603 (2) 102.39 (3) 92.19 (3) 92.05 (1) 104.93 (2) (°) 90 90 91.09(3) 90 90 90 1 (°) -103.3 -95.6 92.0 92.8 -92.6 115.4 2 (°) -178.4 177.5 155.14 159.2 -159.4 -32.0 3 (°) 36.9 -0.9 -10.6 -35.4 8.8 0.0               Graph-set analysis of hydrogen bonds (Etter, 1990; Etter et al., 1990) N8-HO9 C(7) C(7) C(7) C(7) C(7) R22(14) N8-HO10 C(4) C(4) R22(8) R22(8) R22(8) C(4)

Figure 1 Labelling of atoms and torsion angle in piracetam. 1 (torsion angles defined by C2-N1-C6-C7) and 2 (N1-C6-C7-N8) are considered to be the major variable torsion angles, and 3 (C6-C7-N8-H angle closest to 0°) and the torsion angles in the pyrrolidine ring are considered to be the minor variable torsion angles.

During the later stages of this study the authors became aware that form IV of piracetam had been solved (Fabbiani et al., 2005), having been obtained via recrystallization and data collection at high pressure. The experimentalists challenged the authors to test the computational approach by predicting possible structures for form IV, provided only the information that the conformer in the new structure is different from those in the known polymorphs and that the crystal structure was in a common space group with Z' = 1. The challenge was accepted as an excellent opportunity to validate the more subjective aspects of the search process. Six computed structures from the low-energy region (within 5 kJ mol^-1 of the E_tot global minimum) were sent to the experimental team. The six structures were ranked in order of energy, with the lowest-energy most stable structure ranked number 1. This structure proved to be a good approximation to form IV.

2. Methodology

Ab initio calculations using the MP2 level of theory and a 6-31G(d,p) basis set were performed using Gaussian (Frisch et al., 2004) for the conformational optimizations, for the determination of E_intra and to calculate the charge density of the resulting conformers. A distributed multipole analysis (DMA; Stone & Alderton, 1985) for each conformer was obtained from the charge density description using the program GDMA (Stone, 1999). A large number of constrained (con) conformers in which one or more torsion angle was constrained while the rest of the molecule was optimized were considered in this study, along with the gas-phase optimized (opt) conformer and the experimental (exp) conformers taken from the crystal structures with H-atom positions corrected to values expected from neutron diffraction (Allen et al., 1987).

Systematic searches in a number of commonly occurring space groups for low-energy structures were implemented using MOLPAK (Holden et al., 1993) to generate densely packed structures from a pseudo-hard-sphere model of the rigid molecular conformers. In the exploratory searches, eight MOLPAK packing types, corresponding to the P1, P and P2₁/c space groups, were considered. The full searches used 29 MOLPAK packing types, also covering the space groups P2₁, Cc, C2, C2/c, P2₁2₁2, P2₁2₁2₁, Pca2₁, Pna2₁, Pbcn and Pbca. The 50 densest structures in each packing type were lattice-energy minimized using the DMA description of the rigid conformer and the DMAREL algorithm (Willock et al., 1995). The electrostatic contribution to U_latt included all terms in the atom-atom multipole series up to R^-5, with charge-charge, charge-dipole and dipole-dipole terms calculated by Ewald summation. The remaining terms were calculated by direct summation up to a molecule-molecule separation of 15 Å. The empirical repulsion-dispersion potential

was used to model the non-electrostatic contribution, where atom i in molecule 1 is of type and atom k in molecule 2 of type . Parameters for C, N, O, H(-C) (Williams & Cox, 1984; Cox et al., 1981) and polar H atoms, H(-N) (Coombes et al., 1996), were taken from empirically derived potentials.

This model intermolecular potential has been widely used in crystal structure prediction studies and its suitability for piracetam was verified by testing its ability to reproduce the known polymorphs by lattice-energy minimization. All DMAREL lattice-energy minimizations calculate the second-derivative properties, so that any structures that had not reached a true minimum, usually because they were revealed to be transition states, could be discarded. Structures corresponding to lattice-energy minima were sorted by reduced lattice parameters, calculated using PLATON (Spek, 2002), and U_latt to remove exact equivalents.

As the conformer was treated as rigid during each search, two searches using sufficiently similar conformers could result in crystal structures that are effectively equivalent in that they would have led to the same structure if the molecular conformation could have been optimized simultaneously with the lattice energy. Therefore, a means of comparison and classification of low-energy structures from all searches (within a few kJ mol^-1 of the E_tot global minimum) is necessary. The most stable crystal structures were visualized in Mercury (Bruno et al., 2002) for qualitative comparison of packing by consideration of characteristic features in the packing, such as hydrogen-bonding motifs. The similarity of structures with the same motif was quantified using the structure-matching algorithm COMPACK (Chisholm & Motherwell, 2005) for automated pairwise comparison of molecular coordination spheres. The similarity of pairs of crystal structures that matched the non-H intermolecular atom-atom distances in a 15 molecule coordination sphere within (typically) a 20% tolerance was quantified using the r.m.s. deviation of non-H atoms in the optimal overlay of the two 15 molecule clusters.

3. Results

3.1. Exploration of conformational space

3.1.1. Potential energy surface scan

Potential energy surface (PES) scans were carried out about torsion angles ₁ and ₂, both using a 10° step size and considering the full 360° range (see Fig. 2). Each scan was a series of partial, local and gas-phase optimizations with one torsion angle (either ₁ or ₂) constrained and the rest of the molecule unconstrained. Both scans were performed in both directions; discrepancy between two points representing the same results from the dependence of the final conformer on the starting conformer in multiple-minima regions during local optimization. The discrepancy between two minima for a given ₁ typically arises from a difference in ₂ and vice versa, rather than a difference in ring conformation. The clear presence of barriers to conformational change, even in the gas phase, demonstrates one of the challenges associated with CSP for flexible molecules; piracetam has relatively limited flexibility compared with many molecules of similar size, yet the barriers in the PES of ₁ and ₂ (without even considering the other variable torsions) make the location of minima non-trivial.

Figure 2 Ab initio potential energy surface scans of (a) 1 and (b) 2 for piracetam. The scans were carried out in both directions. The discrepancy between the original scans (triangles) and reverse scans (squares) results from the dependence of the final conformer on the starting conformer in multiple-minima regions during the local partial optimization at each point and typically indicates a difference in (a) 2 and (b) 1.

The peaks at 0 and 180° in Fig. 2(a) represent unstable conformations in which the C6-C7 bond partially eclipses the N1-C2 and N1-C5 bonds, respectively. The deep troughs at ₁ = ±90° represent conformers that are stabilized by an intramolecular hydrogen bond between atoms N8 and O9, and the more shallow secondary minimum at ₁ = 90° represents a conformer with a weak intramolecular hydrogen bond between atom O10 and one of the C5 H atoms and no N8O9 interaction. The asymmetry about ±90° for ₁ may be attributed to the asymmetry of the pyrrolidine ring geometry and substituents relative to atom N1. The peak in Fig. 2(b) at ₂ 70° appears to be the result of an eclipse of C6-C7 with N1-C2. At ₂  -90° an intramolecular O9H(-N8) hydrogen bond results in a very stable conformer represented by the wide trough. The trough is likely to be deepened by a simultaneous weak interaction between the O10 atom and a C5 H atom. The discrepancy between structures at ₂ 75° is due to a difference in ₁ (in the more stable of the two minima ₁ = 11°, in the other ₁ = 105°).

₁ and ₂ are considered the major variable torsion angles in piracetam, but all the variations in the molecular geometry will affect the crystal packing to some degree. During the PES scans the pyrrolidine ring varies between the envelope conformation on C4 and the twisted conformation on C4-C5 or C3-C4, and the N1 atom appears to move away from the ring plane when the O10 atom approaches the density of the ring. The NH₂ group tends towards pyramidal geometry to prevent unfavourable close intramolecular contacts (e.g. between protons of N8 and C5) and to facilitate the favourable intramolecular N-HO hydrogen bond, and is more planar in the absence of close contacts.

3.1.2. Comparison with gas-phase conformation

Gas-phase optimizations were started from a number of different conformations to ensure that the global minimum in the MP2 energy was found. The resulting gas-phase optimized conformer has an intramolecular hydrogen bond and is different from the observed conformers in the crystal structures, which are similar in forms II and III but distinct in forms I and IV (see Fig. 3). The side chain is bent in the gas-phase conformer to allow an intramolecular hydrogen bond between the N8 and O9 atoms which is absent in all polymorphs, and in forms I, II and III the N8O9 distance is maximized. The higher level of theory used in the current study means that a discrepancy between the conformation of the side chain in a previously reported (Céolin et al., 1996) semi-empirical AM1 gas-phase optimized molecule (₁ = 74° and ₂ =  298°  =  -62°) and that calculated here (₁ = 80° and ₂ = -81°) is not unexpected. The pyrrolidine ring is twisted (on C4-C5) in the gas-phase optimized conformation, as in forms II and III, yet form IV and both components of disorder in form I have envelope (on C4) conformations. Note that the non-planarity of the ring implies that only the completely inverted molecule is equivalent in terms of energy and packing.

Figure 3 An overlay of the ab initio gas-phase optimized conformer (green) with the observed conformers in form I (major component of disorder only, red), form II (blue, similar to form III) and form IV (yellow). The r.m.s. deviation for all atoms from gas-phase optimized conformer and form I (powder) is 1.96 Å, form I (major component of disorder from single-crystal determination) is 1.91 Å, form II is 1.30 Å, form III (BISMEV01) is 1.29 Å, form III (BISMEV02) is 1.91 Å and form IV is 1.36 Å.

3.2. Reproduction of observed crystal structures as a test of the sensitivity of the lattice-energy minima to the conformation

The good reproductions of the single-crystal structures of piracetam in the CSD by lattice-energy minimization, using the respective experimental conformers (expminexp structures, see Table 2), indicate that the intermolecular potential is adequate, although thermal effects have only been very indirectly included via the use of empirically fitted repulsion-dispersion potentials. The r.m.s. deviations in cell lengths for these reproductions are less than 0.09 Å. The considerable differences between observed conformers and the optimized gas-phase conformer (see Fig. 3) imply that the crystal structures cannot be well modelled using the gas-phase conformer in the lattice-energy minimization, and the resulting expminopt structures were expected to be poor. This was indeed the case, and even the best expminopt reproduction, for BISMEV02, has a 0.85 Å r.m.s. deviation in cell lengths. All other structures have at least three cell parameters in the lattice-energy minimized structures that are shifted by more than 10% from the observed parameters.

Table 2 Reproductions of observed crystal structures of piracetam, using the experimental conformation (expminexp); using the constrained, partially optimized conformation (expmincon) and structures found during the searches with other partially optimized conformations Experimental lattice parameters and, in italics, lattice parameters for calculated structures given as a percentage difference from the experimental values. Form I was found outside the 5 kJ mol-1 range considered, so the `lowest Etot match' (using a 50% tolerance level in COMPACK) outside the 5 kJ mol-1 range is represented rather than the `closest in search' (within 5 kJ mol-1 of the Etot global minimum).               Volume                     per                     molecule Energy# (kJ mol-1)   a (Å) b (Å) c (Å) (°) (°) (°) (Å3) Ulatt Etot F+ Form I (BISMEV03) determined from powder diffraction data Observed 6.747 13.418 8.09 90 99.01 90 180.84 - - - Expminexp 0.04 -0.21 2.43 0 5.92 0 180.96 -109.12 - 93.4 Expmincon 1.33 3.83 -0.85 0 -2.55 0 189.78 -99.82 -86.96 54.1 Lowest Etot match in search§ 2.71 2.06 -11.50 0 -4.21 0 189.41 -99.63 -87.54 -                       Form I (FabbianiI) determined from single-crystal data, only the major component of disorder used in calculations Observed 6.7254 13.2572 8.0529 90 98.603 90 177.48 - - - Expminexp 0.37 -0.57 2.37 0 -0.51 0 181.56 -113.2 - 9.9 Expmincon 2.04 3.54 -0.39 0 -2.36 0 187.77 -99.36 -87.29 73.2 Lowest Etot match in search§ 2.40 0.88 -12.01 0 -4.64 0 189.41 -99.63 -87.54 -                       Form II (BISMEV) Observed 6.403 6.618 8.556 79.85 102.39 91.09 174.25 - - - Expminexp 2.05 -0.45 -1.90 2.72 -0.38 -1.27 174.93 -119.84 - 18.8 Expmincon 0.51 -2.88 7.72 2.87 10.66 0.39 173.32 -104.92 -92.08 256.4 Closest in search 0.25 2.21 -2.14 -20.93 0.15 0.72 175.25 -106.87 -94.06 -                       Form III (BISMEV01) Observed 6.525 6.44 16.463 90 92.19 90 172.82 - - - Expminexp -0.25 0.64 0.81 0 0.89 -0.25 174.80 -118.49 - 3.6 Expmincon -0.99 0.02 2.60 0 6.65 0 173.88 -104.12 -91.97 82.1 Closest in search 1.15 1.68 -4.58 0 0.14 0 175.66 -107.36 -93.65 -                       Form III (BISMEV02) Observed 16.403 6.417 6.504 90 92.05 90 171.04 - - - Expminexp 0.56 -0.08 0.06 0 0.32 0 171.93 -120.00 - 1.2 Expmincon 3.17 -0.31 -0.47 0 4.11 0 174.31 -104.08 -91.80 41.9 Closest in search## -4.96 1.33 0.83 0 -0.01 0 175.66 -107.36 -93.65 -                       Form IV (FabbianiIV) high pressure structure, analysed after prediction work Observed 8.9537 5.4541 13.61 90 104.93 90 160.55 - - - Expminexp 0.04 5.01 2.41 0 1.36 0 171.52 -115.90 - 39.0 Expmincon 1.43 5.87 2.91 0 2.14 0 175.43 -102.36 -91.04 60.5 Closest in search -1.79 -5.36 -2.96 0 -1.95 0 175.48 -103.35 -91.34 - #Eintra and hence Etot cannot be meaningfully calculated for expminexp calculations as the experimentally insignificant changes in bond lengths and angles contribute to the total MP2 energy. +Figure of shame (Filippini & Gavezzotti, 1993) increases with the difference between crystal structures. §Change of cell setting. ##a and c axes swapped.

A series of constrained gas-phase optimizations was performed in which ₁ and ₂ were fixed at experimental values while the remainder of the molecule was relaxed. The resulting rigid conformers were optimally overlaid on the molecules in the experimental crystal structures and lattice-energy minimized (expmincon calculations, see Table 2) to understand the importance of the other unconstrained torsions (₃ and those in the ring) in the crystal packing. These reproductions are typically considerably worse than for the expminexp calculations, indicating that torsions other than ₁ and ₂ are also important in determining the crystal packing, but much better than for the expminopt calculations as the conformation of the side chain (₁ and ₂) is forced to resemble the experimental conformation, confirming that ₁ and ₂ play a major role in determining the packing. The worst reproduction is for form II (0.40 Å r.m.s. deviation in cell lengths), suggesting that torsions other than ₁ and ₂ may play a more important role in the molecular packing in this crystal structure.

The relatively poor expminexp reproduction of form IV, with 0.23 Å r.m.s. deviation in cell lengths, may be explained by the fact that this structure was solved using data collected from a crystal under high pressure. The experimental density is 7% higher than the calculated density, while densities for all other expminexp reproductions agree to within 2.3%. The reproduction of b is poor because of an apparent contraction in this direction in the high-pressure structure (not reproduced in the calculated structure), resulting from the concertina contraction of the zigzag C(4) chain (see Fig. 4) and a corresponding reduction in the size of the void between adjacent pairs of R₂²(14) dimers. The R₂²(14) dimers themselves appear to resist contraction under pressure and they lie in the a direction, which is reproduced more accurately. Discrepancies between the high-pressure structure and the expmincon reproduction follow the same trends.

Figure 4 Overlays of observed high-pressure form IV piracetam (green) with the expminexp reproduction (blue) showing the contraction of the C(4) chains resulting in poor reproduction in the b direction. The R22(14) motif is also shown. H atoms have been omitted for clarity.

It is notable that the structure corresponding to the major component of disorder in FabbianiI is reproduced better than the ordered structure of form I determined from powder data (BISMEV03) in the expminexp calculation. This result may be attributed to the geometry of the amine group in BISMEV03, where the H-atom positions were calculated and not refined, giving a misalignment of NH₂ H atoms with adjacent hydrogen-bond acceptors on neighbouring molecules (N-HO angles 122° and 148°).

The E_tot and U_latt values shown in Table 2 indicate the metastability (and relatively low density) of form I. Forms II and III are close in energy and it is not possible to confidently deduce the order of stability from these calculated U_latt values from either the expminexp or the expmincon calculations. The lattice energy of the newly discovered form IV indicates that it is more thermodynamically stable than form I and less stable than forms II and III. The discrepancies between lattice-energy minimizations for pairs of experimental determinations of the same structure (BISMEV03 and FabbianiI for form I, and BISMEV01 and BISMEV02 for form III) indicate a measure of the sensitivity of the lattice energy to molecular conformation (Beyer et al., 2001); the experimental conformers differ more than the partially optimized conformers and therefore the discrepancies are greater for the expminexp calculations than for the expmincon calculations. The lattice-energy calculations correspond to 0 K and thermal effects are therefore also likely to contribute to the discrepancies for the determinations of form I (one was carried out at 150 K and the other at room temperature).

3.3. Search for low-energy crystal structures

The gas-phase optimized conformer of piracetam seems unlikely to form optimal intermolecular hydrogen bonds in the solid state because of the intramolecular hydrogen bond. The most stable crystal structure found in a search using the gas-phase optimized conformer has E_tot  = U_latt =  -89.6 kJ mol^-1 and has some intermolecular hydrogen bonding. Breaking the intramolecular hydrogen bond in the gas-phase conformer provides the possibility for alternative intermolecular hydrogen bonds to form, in which case a lower U_latt could compensate for a considerable E_intra, and some conformations with very large E_intra were considered. During the early stages of the study, two sets of searches were carried out to probe systematically, but crudely, the crystal packing as a function of ₁ and ₂. Two sets of partial optimizations were carried out; one in which ₁ was constrained and another in which ₂ was constrained (with the remainder of the molecule optimized at each point). The angles were constrained at 20° intervals between 0 and ±160 and 180°. Most of these conformers had no intramolecular hydrogen bond, thus allowing the hydrogen-bond donor and acceptor to be involved in lattice-stabilizing intermolecular hydrogen bonds instead. Each conformer was then used in a search of three common space groups; P1, and P2₁/c (eight MOLPAK packing types). This approach allows the exploration of intermolecular interactions resulting from a range of symmetry elements, and in this case covers the reported polymorphs. The 50 densest structures from each packing type were lattice-energy minimized, resulting in up to 400 structures from each search.

The results of these initial two systematic sets of exploratory searches are represented in Fig. 5; Fig. 5(a) shows E_intra of each conformer, Fig. 5(b) shows the global minimum U_latt from each search, and Fig. 5(c) shows the global minimum E_tot from each search, as a function of ₁ and ₂. Fig. 5 allows quick identification of crude regions in conformational space where there are conformers capable of forming stabilizing intermolecular interactions (such as particular hydrogen-bond motifs), resulting in crystal structures with low E_tot. Comparing the plots in Fig. 5 it is clear that the most stable gas-phase conformers do not form the most stable crystal structures for piracetam. Polymorphs I, II and III are in the region of the low-energy (E_tot) structures to the top right of Fig. 5(c) (₁ 90° and ₂ 170°); low-energy structures elsewhere in the plot (₁ 120° and ₂ -30°, and ₁ 70° and ₂ -160°), therefore correspond to other energetically feasible conformational polymorphs.

Figure 5 (a) Eintra, (b) global minimum Ulatt from the search and (c) global minimum Etot from the search as a function of 1 and 2. Each colour represents a 5 kJ mol-1 range, indicated in the key to the right of the plots; red points indicate the most stable structures, then orange, yellow, green, light blue, dark blue and dark purple, with light purple indicating the least stable structures.

Many other partially optimized conformations (with different numbers of constraints) were used in exploratory searches beyond those represented in Fig. 5. Some conformers were chosen by examining the crystal structures of previous searches and using crystallographic insight to notice that packing motifs and interactions would be improved (or made possible) if the conformer was altered in a specific way (e.g. making a hydrogen-bond donor/acceptor more exposed to allow the formation of a hydrogen bond). Other conformers were chosen simply out of curiosity. To simplify the optimizations and also cover as much of the (₁, ₂) torsion space as possible during the searches, the pyrrolidine ring and ₃ were usually unconstrained. The considerable effect of the orientation of the terminal NH₂ group was investigated in a few cases where ₃ was constrained to, for example, -5° rather than the optimized value of typically -15°. An appreciable effect on E_tot was observed, along with a small increase in E_intra and a considerable decrease in U_latt from an improved crystal packing. However, effectively equivalent crystal structures were generally found (at higher E_tot) for similar conformers with unconstrained ₃; small changes in ₃ are therefore important for the accurate calculation of E_tot but do not give genuinely distinct crystal structures.

Towards the end of the procedure an effort was made to investigate some of the unexplored (₁, ₂) regions, considering some of the blank regions in Fig. 5. These searches produced no better low-energy crystal structures, indicating that the combination of systematic and non-systematic searches had probably generated the majority of low-energy structures that could be determined using this method. In total, approximately 100 searches with different rigid conformers were performed.

It was clear from early in the searching that experimental forms I, II and III had been successfully located, within a well defined region of conformational space, and these conformers were not subject to further searching. The focus of the search was shifted towards different regions of conformational space where other low-energy crystal structures were situated. As a result, six markedly different conformers that produce low-energy crystal packings (crystal structures in the best 5 kJ mol^-1 of E_tot) in the exploratory searches were subject to a full search by considering an additional ten common space groups in case more stable structures could be found with different combinations of symmetry elements. However, only one of these additional space groups produced a new structure within 5 kJ mol^-1 of the global minimum.

3.4. Analysis of low-energy crystal structures

The low-energy crystal structures (within 5 kJ mol^-1 of the global minimum of E_tot = -95.9 kJ mol^-1) are shown in the plot in Fig. 6 and listed as supplementary information. Structures were classified in terms of packing motif and the lowest-energy structure from each packing motif is given in Table 3. The penultimate column in Table 3 is the packing motif label; observed forms of piracetam II, III and IV are labelled as such, while other motifs (1-10) are hypothetical structures that are, as yet, unreported experimentally. The structures listed in Table 3 are provided (in CIF format) as supplementary information.¹ The classification of structures is based on the 20% distance tolerance level in COMPACK (the r.m.s. deviation between structures with the same motif is typically around 0.3 Å). If this tolerance is increased there are fewer distinct motifs in terms of a 15 molecule coordination sphere. Motif 2 is the same as II at 50% tolerance (r.m.s. deviation is typically 1.5 Å), motif 4 is the same at 80% tolerance (r.m.s deviation typically > 1.5 Å) and motif 5 is the same again at 100% tolerance (r.m.s. deviation of the most stable motif 5 structure and experimental form II is 2.07 Å). Similarly, motif 1 is the same as motif 10 at 80% tolerance (r.m.s. deviation typically 1.2 Å). Other motifs remain distinct up to 100% tolerance. This table illustrates the challenge of deciding how similar two computed structures must be to be likely to correspond to the same structure experimentally, either through the realistic optimization of all the atomic positions producing the same minimum or through easy transformation (perhaps simply by thermal motion) between closely related minima. Perhaps the experimental observation of motif 5 is unlikely because of a facile transformation to form II, for example. Certainly in the case of piracetam the observed polymorphs are distinct up to 100% tolerance.

Table 3 Crystal structures representing the distinct packing motifs in the lowest 5 kJ mol-1 of Etot in the computational search Two crystal structures are listed for each of the observed packing motifs (II, III and IV) with r.m.s. deviation (r.m.s.d., 15 molecule coordination sphere) from best experimental structure; the first is the lowest Etot structure of that type and the second is the closest (in terms of r.m.s. deviation) to the observed crystal structure. The lowest Etot representative from each of the hypothetical structures (1-10) is also listed.             Volume                   per       Conformer     Energy (kJ mol-1) molecule Space Packing r.m.s.d. label 1 (°) 2 (°) Eintra Ulatt Etot (Å3) group motif (Å) 3a 88 160 13.05 -108.95 -95.90 173.93 P21/c III 0.43 3a 88 160 13.05 -108.90 -95.84 173.21 P21/n 1   3a 88 160 13.05 -108.16 -95.11 172.32 2   3c 90 160 13.30 -108.23 -94.93 175.49 II 0.27 3d 90 160 12.81 -106.87 -94.06 175.25 II 0.26 10d 100 160 13.71 -107.36 -93.65 175.66 P21/c III 0.22 10b 120 -34 12.55 -106.10 -93.56 174.6 P21/c IV 0.28 10g 80 160 13.57 -106.73 -93.16 175.38 P21/c 3   1d 69 -160 10.81 -103.16 -92.35 176.03 4   10b 120 -34 12.55 -104.90 -92.35 179.3 5   1d 69 -160 10.81 -102.75 -91.94 179.07 P21/c 6   3c 90 160 13.30 -105.21 -91.91 178.81 P21/c 7   6f 120 -36 12.01 -103.35 -91.34 175.48 P21/c IV 0.27 9o 120 160 15.77 -106.89 -91.12 174.39 P21/n 8   10b 120 -34 12.55 -103.66 -91.11 182.41 C2/c 9   1a 68 -150 11.18 -102.27 -91.09 179.02 P21/c 10

Figure 6 Crystal structures within 5 kJ mol-1 of the global minimum in terms of Etot. Each symbol represents a different structure type as defined in the key. Black points are calculated hypothetical structure types (1-10) and red points are calculated structures types that have been observed experimentally (red squares, diamonds and triangles represent type II, III and IV structures, respectively). Results of expmincon calculations (solid red points, see Table 2) do not correspond exactly to structures found in the searches (hollow red points) because 1 and 2 were constrained to the exact experimental values in expmincon calculations. Form I calculated structures are outside the energy range in this plot.

It is clear from Fig. 6 that there are a number of crystal structures (with a variety of conformers) that are energetically competitive with the observed forms. C(7) and R₂²(14) motifs involving the N8/O9 donor/acceptor pair and C(4) and R₂²(8) motifs with the N8/O10 donor/acceptor pair are the hydrogen-bond motifs found in the low-energy structures. The N8 atom forms bifurcated hydrogen bonds, so Z' = 1 structures of piracetam are expected to be stabilized by a combination of two hydrogen-bond motifs. The resulting four pairwise combinations of hydrogen-bond motifs are energetically competitive and all occur in the 5 kJ mol^-1 range considered here, with the R₂²(8) motif appearing most frequently (in 11 out of the 13 structures). Of the structures with the R₂²(8) motif, structures 4, 5, 6 and 9 also exhibit R₂²(14) dimers, while structures 1, 2, 3, 7, 10, II and III also exhibit C(7) chains. The C(7) chains stack differently with little variation in energy. Consecutive molecules in the linear C(7) chains in structures 1, 2, 3, 7, II and III are related by translation, and all conformers in these structures have ₂ = 160° (with one exception for a type 1 structure just within the 5 kJ mol^-1 range at E_tot = -90.91 kJ mol^-1, where ₂ = 140°) and ₁ within a 20° range. In structures 8 and 10 consecutive molecules in the C(7) chains are related by a glide plane and a screw axis, respectively, resulting in zigzag chains. Structure 8 also has elongated C(4) contacts rather than R₂²(8) dimers, like form IV, which is stabilized by C(4) and R₂²(14) interactions. Clearly the possibility of this range of complex hydrogen-bonding motifs with similar energies shows that conformational polymorphism is thermodynamically feasible. However, there is no clear association between motifs and conformers in the observed polymorphs to suggest a kinetic factor in the early stages of molecular association favouring the metastable polymorphs. It is worth noting that the use of high pressure was necessary to provide the conditions for the conformational polymorph form IV of piracetam to grow.

Forms II and III are reproduced frequently within the low-energy structures. The quality of reproduction of forms II and III in the search is good (Table 3); the best r.m.s. deviations with experimental forms II (BISMEV) and III (BISMEV02) are 0.26 and 0.22 Å, respectively. The reproduction of cell parameters is shown in Table 2. The prediction that form III at the E_tot global minimum, and form II close in energy, are the most stable forms at 0 K is a success for the accuracy of the modelling. Recent work (Jagielska et al., 2004) to derive parameters for a non-bonded potential found piracetam form III (BISMEV02) as rank 9, 4.5 kJ mol^-1 less stable than the global minimum, although many other known amide structures were found as global minima in these crystal structure prediction studies using the rigid experimental conformer. Metastable form I was found in the searches at energies outside the 5 kJ mol^-1 range, as expected from expmincon calculations. The most stable occurrence is at E_tot = -89.1 kJ mol^-1 (nearly 7 kJ mol^-1 less stable than the global minimum); however, this is a poor reproduction with an r.m.s. deviation of 2.06 Å from the observed form (using the major component of disorder for the experimental structure, and ignoring H atoms). A better (but less stable) reproduction occurs at E_tot = -87.5 kJ mol^-1, with an r.m.s. deviation of 0.48 Å (see Table 2).

3.5. Results of blind prediction of form IV

Form IV was determined (Fabbiani et al., 2005) during the later stages of the searching and was known to the authors to be a Z' = 1 conformational polymorph of piracetam. When the searching was complete, the seven packing motifs with markedly different conformers to those previously known, in the lowest 5 kJ mol^-1 of E_tot, were possibilities for the new form IV. Only one structure in the 5 kJ mol^-1 range, representing packing motif 9, could be disregarded because of its low density (182.4 Å³ per molecule compared with 179.3 Å³ per molecule for the next least dense structure in the 5 kJ mol^-1 range); it was considered unlikely to form under high pressure. The remaining six packing motifs were all considered reasonable predictions, so the lowest-energy crystal structure from each was suggested to the experimental team who had independently solved the new structure. These six crystal structures are included in Table 3 and belong to packing motifs IV, 4, 5, 6, 8 and 10 (in order of decreasing stability). The most stable of these crystal structures (packing motif IV), E_tot = -93.56 kJ mol^-1, was confirmed by the experimentalists to correspond to the newly determined form IV (note that the label `IV' was assigned after this confirmation). The r.m.s. deviation between observed and predicted form IV is 0.05 Å for the conformer and 0.28 Å for a coordination sphere of 15 molecules (both ignoring H atoms) and a unit cell overlay is shown in Fig. 7. This study therefore represents a successful `blind' crystal structure prediction for a flexible molecule in which the conformer is markedly different from both the gas-phase optimized conformer and observed conformers in other experimental polymorphs.

Figure 7 Observed unit-cell contents of form IV (green) with the predicted structure (blue) overlaid.

4. Discussion

An efficient approach to CSP using a partial search of conformational space has been used to implement a computational investigation of polymorphism in piracetam. A search involving the ab initio optimization of both conformer and crystal structure is certainly unrealistic in this case (in terms of time and computing resources) where the search space is so vast. Here, rigid conformers with ab initio calculated intramolecular energies were used as rigid probes in the search for stable crystal structures, which were then lattice-energy minimized. It is necessary to achieve an appropriate balance between the systematic investigation of conformers and refinements of promising conformers based on crystallographic insight and to realize that some conformation alterations improve the calculated stability for a crystal structure but do not lead to the computation of novel crystal structures. The coverage of conformational space, the choice of step size and the extent of each search required to investigate a molecular system sufficiently will vary and could be estimated from the energy profiles of the variable torsion angles. The consideration of a large number of molecular geometries was necessary in the search for hypothetical piracetam structures because of the need to explore the different combinations of competing hydrogen-bond motifs. The range of conformations that need to be considered would be smaller for many molecules of similar size and conformational flexibility. For example, if the exploratory searches showed that the molecule had just one possible hydrogen-bond motif, so that low-energy structures were only accessible to a small range of conformations, then only this range of conformers would need investigation in more detail, with the possibility of using more accurate energies. On the other hand, the approach developed in this study is most appropriate for molecules capable of forming hydrogen bonds, and more emphasis would be needed on the systematic component of the search for molecules without such specific interactions.

The style of searching demonstrated here is aided greatly by the availability of grid compute facilities (Butchart et al., 2003). The parallel distribution of MOLPAK and DMAREL jobs over a large number of nodes means that the entire search is completed in a much shorter time (it is realistic to expect an exploratory search to run in under an hour); thus, not only is it possible to consider more conformers, but also, importantly, the choice of rigid conformer can be based on results from previous searches so that the search may be driven `interactively'. Studies of piracetam were conducted using both the traditional approach, in which jobs are run sequentially on a Linux cluster, and new grid compute facilities, in which the jobs are invoked via web services and run in parallel. Even during the testing phase of the grid compute facilities, the benefit they provide for studies such as this was clear, and it is envisaged that future studies will benefit further.

The newly determined form IV of piracetam has been successfully predicted `blind', to a considerable degree of accuracy, despite differences between the form IV conformation and the gas-phase optimized conformation and other observed conformations. There are a number of hypothetical structures that are competitive in terms of energy with the observed polymorphs and therefore could be undiscovered polymorphs, given the lack of understanding of the kinetic factors involved in the crystallization of piracetam. The range of conformations that allow stable crystal packing motifs for this molecule is evident in the low-energy calculated structures. The marked difference between gas-phase optimized and observed conformers demonstrates the considerable thermodynamic benefit for piracetam in terms of intermolecular lattice-stabilizing interactions to be achieved with a suboptimal gas-phase conformer. Extension to other flexible molecules and salts (two fragments in the asymmetric unit) will determine whether this approach is feasible more generally.

5. Conclusions

Flexible molecule crystal structure prediction requires the calculation of close-packed crystal structures where the molecular conformation facilitates stabilizing intermolecular interactions. The successful blind crystal structure prediction of piracetam suggests that the approach described here, which has a number of benefits over alternative methods, is promising for flexible molecules with more than one variable torsion angle and a range of competing hydrogen-bond motifs.