Amide Pyramidalization in Carbamazepine - American Chemical Society

Amide Pyramidalization in Carbamazepine: A Flexibility Problem in Crystal Structure Prediction? Cabeza,†

Aurora J. Cruz William Jones*,†

Graeme M.

Day,†

W. D. Samuel

Motherwell,‡

and

CRYSTAL GROWTH & DESIGN 2006 VOL. 6, NO. 8 1858-1866

The Pfizer Institute for Pharmaceutical Materials Science, Department of Chemistry, UniVersity of Cambridge, Lensfield Road, Cambridge CB2 1EW, U.K., and Cambridge Crystallographic Data Centre, 12 Union Road, Cambridge CB2 1EZ, U.K. ReceiVed March 30, 2006; ReVised Manuscript ReceiVed May 22, 2006

ABSTRACT: Carbamazepine is known to exist in various polymorphic forms. Here we report on crystal structure prediction calculations for carbamazepine in an attempt to examine the predictability and relative stability of the various polymorphs. Hypothetical crystal structures generated in 10 of the most common space groups were compared to the observed polymorphs. Particular attention has been given to the influence of amide pyramidalization on the relative energies of the predicted structures. While the actual generation of possible crystal structures was found to be independent of the degree of deformation of the amide group, their final ranking in energy was greatly affected by pyramidalization of the amide nitrogen. This effect was examined in detail through systematic variation of the NH2 geometry for each of the low-energy crystal structures; different amide geometries were favored in the various low-energy crystal structures. The results demonstrate that energetically feasible deformation of the amide group may occur in order to optimize hydrogen-bond interactions, and we conclude that neglect of amide pyramidalization introduces significant errors in crystal structure prediction for carbamazepine and similar molecules. Introduction Polymorphism, the ability of a substance to crystallize in a variety of different crystal structures, is of considerable interest. For a drug molecule, polymorphism is recognized as an undesirable feature, since the various polymorphs may display quite different physical properties.1-3 The ability to predict the relative thermodynamic stability and likely properties of known and latent polymorphs is, therefore, clearly of importance to various applications of molecular solids, including pharmaceuticals. Such predictions, however, remain a formidable challenge, with many difficulties to be overcome before reliable predictions can be applied to typical pharmaceutical systems.4-7 One particular problem is in the development of strategies and methodologies for dealing with conformationally flexible molecules.8,9 Carbamazepine (CBZ; Figure 1), an anticonvulsant drug used in the treatment of epilepsy and trigeminal neuralgia, has, to date, been observed to crystallize in four polymorphic forms (forms I-IV, Figure 2),10-13 with the most recently observed form being characterized only 2 years ago (despite 30 years of study of this system). Concurrent with the work described here, automated parallel crystallization methods have been used in conjunction with crystal structure prediction, to explore the polymorphism of carbamazepine under a wide range of crystallization conditions.14 No new polymorphs were observed after 594 individual crystallizations. Additionally, although calculations suggested that structures based on a chain hydrogen bond motif were energetically favorable, the experimental observations showed that, during crystallization, formation of the R22(8) carboxamide dimersthe motif observed in all four polymorphsswas preferred.14 Here we report a more detailed computational investigation of carbamazepine, with particular emphasis given to the influence of amide pyramidalization on the relative energies of the real and hypothetical crystal structures. Although amide * To whom correspondence should be addressed. E-mail: [email protected]. † University of Cambridge. ‡ Cambridge Crystallographic Data Centre.

Figure 1. Molecular structure of carbamazepine.

Figure 2. Unit cells of the experimentally known carbamazepine polymorphs.

groups are generally planar, this planarity can be distorted through C-N bond rotation (twisted amides) or nitrogen pyramidalization (pyramidal amides). Nitrogen pyramidalization is energetically more feasible than C-N bond rotation but, in general, both processes are very closely related and neither occurs without some contribution of the other. This is a wellknown problem in the computer simulation of peptides.15 Although certain force fields attempt to address this matter,12,16 the problem is complicated by the fact that a change in nitrogen pyramidalization (and hence hybridization) will be accompanied

10.1021/cg0601756 CCC: $33.50 © 2006 American Chemical Society Published on Web 06/27/2006

Amide Pyramidalization in Carbamazepine

Crystal Growth & Design, Vol. 6, No. 8, 2006 1859

by large changes in the partial atomic charges typically used to model intermolecular electrostatic interactions. As a result, conformational dependence of the electrostatic model must be considered in order to correctly describe the molecular electrostatic potential. The sensitivity of crystal structure modeling and prediction to small changes of geometry of an amide group in an otherwise rigid molecule has not hitherto been addressed, due in part to a lack of accurate experimental data. Hydrogen atoms, often the most important atoms involved in intermolecular interactions in molecular solids, are the most poorly resolved by X-ray diffraction, the most routinely used technique for crystal structure determination. We report here on our studies concerning the predictability of three (forms II-IV) of the four polymorphs of carbamazepine. We exclude form I, which has four molecules in the asymmetric unit and is beyond current methods in crystal structure prediction, and limit our studies to crystal structures with Z′ ) 1. We compare two descriptions of the electrostatic intermolecular interactions: (i) molecular electrostatic potential (ESP) derived atomic point charges17,18 and (ii) atomic multipoles obtained from a distributed multipole analysis.19 In addition, the influence of the inversion of the amide nitrogen on the final energies of the known and hypothetical polymorphs is examined. Since adequate force fields with suitable terms for the description of the NH2 pyramidalization, together with an accurate description of intermolecular interactions, are not available, we are limited to treating the molecules as rigid bodies in a series of calculations using systematically varied molecular geometries, derived from isolated molecule first principles calculations.

density obtained using the CADPAC code.26 All charge-charge, charge-dipole, and dipole-dipole interactions were evaluated using Ewald summations, while higher order electrostatic terms (up to R-5) were summed to a 15 Å cutoff on the basis of whole molecules. Procedure for Crystal Structure Prediction. Hypothetical crystal structures were generated in the nine most common space groups (P21/c, P1h, P212121, P21, C2/c, Pbca, Pnma, Pna21, and Pbcn) using the simulated annealing algorithm of Karfunkel and Gdanitz,27-29 as implemented in the Accelrys Polymorph Predictor (PP) module of the Cerius2 software suite.30 Adjustable parameters controlling the simulated annealing were taken from previous studies.31 In addition, we searched space group R3h, as one of the known polymorphs exhibits that crystal symmetry. A blind crystal structure prediction would, of course, need to consider a fuller set of space groups; here, we limit ourselves to the nine most common, plus R3h, to save on computational cost, while searching the space groups occupied by the known carbamazepine polymorphs and the most likely alternate crystal packings. Consecutive independent searches from random starting points (up to 3 or 4) were performed within each space group until additional structures were no longer generated. These initial searches were carried out with one molecule in the asymmetric unit (Z′ ) 1) and were followed by lattice energy minimization (W99 + ESP charges). Hypothetical structures were then clustered (using the algorithm implemented in the Cerius2 software29) to remove identical structures. Lattice energy minimization of the predicted structures was also performed with the atomic multipole electrostatic model and the same W99 model potential using the program DMAREL.32 Final Clustering. Predicted structures were finally clustered using the COMPACK33 algorithm with a tolerance of 0.15 Å in atomic separations (the COMPACK algorithm compares structures on the basis of atomic separations in a coordination sphere of 15 moleculessa central molecule and its 14 neighbors). The same algorithm was used for comparison of predicted structures with the experimental structures.

Computational Methods

Exploring the Carbamazepine Molecular Energy Surface and Choosing a Starting Conformation. Prior to searching for crystal structures, it is necessary to identify possible lowenergy molecular conformations. Gas-phase geometry optimizations of the isolated carbamazepine molecule were therefore performed. Two molecular conformations were identified as minima in the potential energy surface (Figure 3). The global minimum (denoted as Hin) was about 2 kJ/mol more stable than a second, local minimum (designated Hout). The only significant difference between the two geometries is the direction of the nitrogen pyramidalization and, as a consequence, the orientation of the hydrogen atoms of the amide group. Using the OCN molecular plane as reference, the amide hydrogen atoms may be described as pointing toward the inner (Hin) or the outer (Hout) part of the molecule: see Figure 3. The lowest energy interconversion pathway between these two minima was studied using Dmol3 and the relative conformational energies obtained as a function of the selected HNCN dihedral angle (Figure 4); calculations were performed by constraining the HNCN dihedral angle at several values from -60° to +50°, optimizing the remainder of the molecule at each point. Each value of the HNCN dihedral angle corresponds to a different degree of amide pyramidalization. Negative values result in the amide hydrogens pointing below the plane (Figure 3) and positive values above the plane: 0° corresponds to the planar geometry. As expected, the energy barrier between molecular minima is sufficiently small (2.06 kJ/mol) to allow interconversion in the gas phase, although in the crystal, hydrogen bonding may limit the mobility of the hydrogen atoms. Rotation about the C-NH2 bond was also monitored but was found to be significantly more expensive in terms of energy (the calculated rotation barrier was almost 46 kJ/mol). Nonbonded intramolecular interactions of the amide group (CONH2) with the rest of the molecule were also evaluated for

Molecular Model. All molecular models were treated as rigid during the crystal structure search and lattice energy minimizations. Initial molecular structures were obtained using density functional theory (DFT) geometry optimizations of the isolated molecule using the code Dmol3 (as implemented in the Accelrys package Materials Studio20) and the PW91 functional21 along with a double numerical polarized (DNP) basis set.22 In studying the effect of nitrogen pyramidalization, conformational energies were evaluated at the same level of theory by optimizing molecular structures with HNCN torsion angles constrained at a series of values. Potential Model. The repulsion-dispersion contributions to the intermolecular potential were evaluated as Nmol

Urep-disp )

∑( ∑

Uikvdw) )

M,N i ∈ M < k ∈ N Nmol

(

∑ ∑

M,N i ∈ M < k ∈ N

(

Aικe-BικRik -

Cικ

))

Rik6

where atoms i and k in molecules M and N are of types ι and κ, respectively, and the parameters Aικ, Bικ, and Cικ are characteristic of the atom types. The Williams W99 empirically derived atom-atom potential parameters23-25 were used, and interactions were summed to a 15 Å cutoff. In this potential, the center of interaction for all hydrogen atoms was shifted by 0.1 Å along the X-H bond (X ) C, N, O) toward the heavy atom and, therefore, the X-H distances were shortened by 0.1 Å after DFT optimization. Electrostatic Models. Electrostatic contributions to the intermolecular potential were calculated using two models for the molecular charge distribution: atom-centered charges and atomic multipoles. Atom-centered point charges were derived to reproduce the molecular electrostatic potential (ESP charges)17,18 using the fitting procedure implemented in Dmol3. Atomic multipoles up to l ) 4 (charge, dipole, quadrupole, octupole, and hexadecapole) were obtained by performing a distributed multipole analysis19 of the B3P91/6-31G(d,p) electron

Results

1860 Crystal Growth & Design, Vol. 6, No. 8, 2006

Cruz Cabeza et al.

Figure 3. Gas-phase carbamazepine molecular geometry minima and planar geometry and their relative energies (∆EbondedQM) referenced to the most stable conformation (Dmol3, pw91/dnp).

Figure 4. Hin-Hout gas-phase interconversion pathway (Dmol3, pw91/ dnp).

the conformations along the interconversion pathway (using the W99 parameters for the repulsion-dispersion terms and the atomic point charge model for the electrostatic component). Subtracting this nonbonded contribution from the DFT values calculated gave a parabola with minimum at 0°. These calculations suggest that the isolated amide group prefers a planar conformation and that the amide deformation from planarity in carbamazepine is a consequence of the nonbonded intramolecular interactions of the hydrogen atoms with the rest of the molecule. These interactions are of the same type and magnitude as those between molecules in the crystal and, therefore, crystalpacking forces will have a significant influence on the degree of deformation of the amide group. As a result, it is important to investigate how the choice of a starting molecular model influences the outcome of the various steps involved in crystal structure prediction: i.e., the generation of all closely packed structures and the evaluation of their relative stabilities. Although the conformation adopted by the molecule in the crystal does not necessarily correspond to an energy minimum for the isolated molecule, we first performed calculations with the Hin

and Hout minima to study how the choice of conformation affects the generation of crystal structures and their relative energies. Effect of Molecular Conformation on Crystal Structure Generation and Ranking. The Monte Carlo simulated annealing procedure was used to generate crystal structures using both the Hin and Hout conformations and the relative lattice energies (∆Enon-bonded ) Enon-bonded,i - Enon-bonded,1, referenced to the most stable crystal structure, ranked i ) 1), and densities of the structures were analyzed after lattice energy minimization with the atomic point charge and the atomic multipole models (Figure 5). As stated earlier, prediction of the Z′ ) 4 form I structure was not attempted. Instead, starting from the experimental crystal structure, form I was lattice energy minimized after replacing the X-ray-determined molecular geometries by the Hin and Hout molecular models (assuming that all four independent molecules adopt the same geometry). We found that (in the low-energy regions of packing space, where the simulated annealing sampling is complete) the same set of hypothetical crystal structures was generated with both molecular geometries; i.e., there is a one-to-one correspondence of structures generated with Hin with those generated with Hout. The three Z′ ) 1 experimental polymorphs (circled structures) were generated in both searches; the conformational differences between the molecular models were not sufficiently important to affect the generation of crystal structures. Nevertheless, the position of the NH2 hydrogen atoms clearly affected the calculated lattice energies (Figure 5) and, therefore, the ranking of the structures. This result is unsurprising: these polar hydrogen atoms contribute very little to the overall molecular shape, which guides the search for close-packed crystal structures, but play a major role in the hydrogen bond interactions that are important contributors to the overall lattice energies and, hence, the relative stabilities of the crystal structures. In the Hout set of predictions, the first- and second-ranked hypothetical crystals with both electrostatic models correspond to the experimental monoclinic structures (forms III and IV), whereas with Hin an experimentally unobserved structure was found to be the most stable. Furthermore, it is to be noted that the order of stability of forms III and IV for the two models is reversed, depending on the conformation used. Comparison of the total energy (Eitotal ) Enon-bonded,i + ∆EbondedQM) of the lowest energy predicted crystal structures (those found within 3 kJ/mol of the global energy minimum for each molecular model) shows that the choice of amide conformation has a different influence on each crystal structure



Figure 5. Relative energies and densities of predicted crystal structures of carbamazepine obtained with Hin (diamonds, above) and Hout (squares, below) molecular conformations, after minimization using the atomic point charge model (left) and the atomic multipole model (right). Circled structures correspond to those matching the experimentally observed polymorphs. Form I is lattice energy minimized from the X-ray determined structure. Table 1. Two Structural Angles That Define the Conformations Used to Reminimize the Computer-Generated Crystal Structuresa conformation

HNCN (deg)

H′CNH (deg)

1 2 (Hin) 3 4 (planar geometry) 5 (Hout) 6

-30.0 -24.3 -15.0 0.0 14.6 30.0

-136.9 -141.8 -150.7 180.0 153.4 138.9

a H is the amide hydrogen that points toward the aromatic ring system, whereas H′ points away from it. HNCN is a proper torsion, whereas H′CNH is improper (out of plane).

Figure 6. Total energy shift of the same crystal structures modeled using Hin and Hout molecular conformations (W99 + multipole model). Solid lines join crystal structures containing dimers, whereas structures joined by dashed lines contain chain motifs.

and can change the total relative energies by more than 8 kJ/ mol (Figure 6). All of the structures in this energy window contain either hydrogen-bonded dimers or chain motifs. The lines indicate matching crystals with Hin and Hout molecular

structuresssolid lines join those crystals with hydrogen-bonded dimers, and dashed lines join those structures containing chains. These shifts in relative energies seem to indicate that amide pyramidalization has a greater effect on the energies of structures involving hydrogen-bonded chains, as compared to its effect on dimers. Influence of Amide Pyramidalization on the Relative Total Energy and Ranking of the Predicted Carbamazepine Crystals: Refining the Calculations. The dependence of the lattice energy on the NH2 geometry suggests the need to refine the calculations in order to consider the possible deformation of the amide group from the two gas-phase minima within the crystal. A selection of molecular conformations along the HinHout interconversion pathway (Figure 4) was therefore considered, and the previously generated set of crystal structures was reminimized after successive incorporation of each molecular conformation. Six molecular conformations (Table 1) were selected for these calculations.


Figure 7. Influence of the molecular amide pyramidalization on the relative total energy of the predicted crystal structures using the atomic point charge model. The thickest lines correspond to the structures matching the four observed polymorphs.

For each crystal structure generated by the simulated annealing searches with the Hin and Hout molecular models, the original molecular model (either Hin or Hout) was replaced by each of those listed in Table 1, aligning all the atoms other than the amide hydrogens. Each predicted crystal structure was then reminimized using the W99 + atomic point charge model. The relative total energy was calculated as a sum of the intermolecular and conformational energies (intermolecular energies calculated with the W99 + atomic point charge model and conformational energies from DFT):

∆Etotal ) Eitotal - E1total ) (Enon-bonded,i + ∆EbondedQM) - (Enon-bonded,1 + ∆EbondedQM) and the total energy for each crystal structure was calculated as a function of the HNCN dihedral angle (Figure 7). In some cases (e.g. forms I and II), when the lowest energy point was located at the limit of the chosen conformational range, an additional conformation was considered (HNCN ) 45.0°, H′CNH ) 128.4°). Each of the chosen conformations was found to provide a minimum energy in one or more of the predicted structures, with many of the low-energy crystals favoring a molecular geometry away from the gas-phase minima. The energies of some crystal structures, such as the lowest energy hydrogenbonded chain structure (Figure 7, triangles), were greatly affected by changes in conformation, while others, such as form III, show a greater freedom of molecular deformation with small

Cruz Cabeza et al.

accompanying energetic cost. The importance of molecular flexibility in the ordering of the possible low-energy structures is clear; no single molecular geometry gives the correct energetic ordering of the different crystal structures. Allowing molecular deformations is therefore essential to obtain the correct crystal stability. Similar calculations were also performed with the atomic multipole model for a subset of these crystal structures (see the Supporting Information), and we found that the favored amide geometry for each crystal structure did not depend on the electrostatic modelswhen multipoles were used, the crystal energy curves were shifted in energy but with the same shape and minimum energy dihedral angle as with the atomic charge model (Figure 7). Because the minimum energy conformation was unaffected by the electrostatic model, the dihedral angle scan did not have to be repeated with atomic multipoles for the entire set of crystal structures, saving the large computational expense associated with such calculations. Instead, atomic multipoles were only used in the final step in our calculations, where each crystal structure was reminimized with the W99 + atomic multipole model (Figure 8), keeping the molecular geometry for each crystal structure fixed at the minimum energy conformation found from the W99 + atomic point charge dihedral angle scans (Figure 7). The atomic multipole description of the electrostatics improves the ranking over the atomic point charge model for two of the observed polymorphs of carbamazepine: forms III and IV (the two monoclinic structures) are calculated to be two of the three most stable structures, together with a third form, which has not been experimentally observed (Figure 8). Forms I and II (the triclinic and trigonal polymorphs) were not so well ranked on total energy, being 6.2 and 8.3 kJ/mol above the global minimum; if they had not been observed experimentally, the calculations described here would not have predicted either of them as likely polymorphs. Discussion Predicted Stability of the Observed Polymorphs. We have identified several factors that are important for the modeling and prediction of the crystal structures of carbamazepine. As demonstrated in earlier work,34 the form of the electrostatic model has a significant impact on the lattice energy ranking of computer-generated hypothetical crystal structures, with the theoretically better founded atomic multipole model typically yielding more reliable results (the experimentally observed structures closer to the global minimum in lattice energy of the computer generated crystal structures). Changes to the geometry of the amide group also have an important effect on the relative stabilities of the observed and low-energy hypothetical structures. Our final model finds two of the known polymorphs (III and IV) among the three lowest energy predictions, while the other two known polymorphs (I and II) are ranked higher in energy than many alternative hypothetical polymorphs, at more than 6 kJ/mol above the global minimum. All calculations presented here are static and temperature-free, and it should be noted that there is a 10% variation in densities (1.19-1.32 g/cm3) among the set of low-energy structures; this could lead to important differences in thermal contributions to the relative free energies at temperatures above 0 K. Lattice dynamics calculations might, therefore, significantly alter the energy ranking of the real and hypothetical crystals. Form I (P1h) is obtained from the melt and shows many similarities to form II (R3h). It is known that form II transforms into form I upon heating. Our calculations place form I at 1.73



Figure 8. Final carbamazepine crystal structure predictions with atomic point charge (left) and atomic multipole (right) models, after the dihedral angle scan to locate the preferred amide conformation in each crystal structure. Structures are labeled by their hydrogen bond motifs: dimers (circles), chains (squares), and structures with no conventional hydrogen bonds (triangles). Full structures of the lowest energy crystal structures are available in the Supporting Information.

kJ/mol lower energy than form II, an energy difference which is within the range of the experimentally observed exothermic form II f I transition enthalpies (measured values range from -1.46 to -2.34 kJ/mol).13,35 Nevertheless, we calculate form I to be less stable than form III (P21/c) by 6.2 kJ/mol, which is much higher than the experimentally observed transition enthalpy (between 1.30 and 3.60 kJ/mol).13,36-40 The form II f form I transition enthalpy can be measured directly with DSC, whereas ∆H for form III f form I is always obtained indirectly, with little agreement among the values reported in the literature. In comparing the experimental and theoretical data, we must recognize the possible errors that affect both the measurements and calculations. Experimental accuracy might be poor if multiple processes are involved during the DSC measurements, while the accuracy of the calculations is limited by the use of empirical atom-atom repulsion-dispersion potentials, the static energy minimization approach, and the uncertainties in the molecular geometry. In this study, molecular flexibility due to amide pyramidalization has been considered carefully but it should be noted that other sources of molecular flexibility might be possible,41 including pyramidal distortion of the ring nitrogen, breathing of the azepine ring, and tilting and rotation of the carboxamide group. Rotation about the N-CONH2 bond in the isolated carbamazepine molecule was investigated but was judged to be less important than NH2 pyramidalization and so was not considered in the crystal energy calculations. It is also important to point out that the gas-phase optimized molecular structure agreed better with the experimental molecular geometries found in form III (CBMZPN02) and form IV (CBMZPN12) than those in form II (CBMZPN03) or form I (CBMZPN11), where the angle between rings in the azepine ring system is distorted. By ignoring this type of molecular flexibility, we will have overestimated the relative energies of forms I and II, which must be part of the reason for their poor ranking in lattice energy. Nevertheless, a qualitative agreement in the relative stability of the solid modifications should be expected.42 Our final calculations give the following order of stability at 0 K: form III (P21/c) ≈ form IV (C2/c) > form I (P1h) > form II (R3h), in good agreement with experimental studies (form III > form I > form II). The only polymorph in disagreement is form IV, which is reported13 to lie between

forms I and II. Only one experimental investigation of the stability of this polymorph, however, has been reported. Molecular Conformations in the Crystals. Different degrees of amide pyramidalization were found in the final energy minimized structures for forms II-IV (Table 2). The final minimum for form I has the same molecular geometry as form II, though independent variation of the geometry of the four molecules within the asymmetric unit was not allowed in the calculations. We note that the reported HNCN torsion angles in form III (+2.5, -4.0, and -6.3° in the three determinations) are all close to the planar geometry that we calculate as the energy minimum, while the amide geometries in the X-raydetermined form I structure (HNCN angles of +7.0, +11.3, +15.8, and +21.7° for the four molecules in the asymmetric unit) are all further from planar; in this case, our calculations also gave a significantly nonplanar geometry of 30°. Experiment (HNCN ) +3°) and calculation (15°) are in less good agreement for the C2/c polymorph, and no hydrogen positions were experimentally determined for the R3h structure. However, because of the limitations of locating H atoms accurately from X-ray data, such comparisons are not a reliable assessment of the computationally derived amide geometries; neutron diffraction studies would be required to assess the computational results. Although all four polymorphic forms involve hydrogenbonded dimers, each exhibits small structural differences in the orientation of the two molecules (Figure 9)sthe NCO planes of the hydrogen-bonded molecules are not coplanar in all four polymorphs. In form III, the NCO planes are aligned, leading to optimal hydrogen bonding with the planar NH2 geometry. However, planar amide groups in form IV give a nonlinear N-H- - -OdC dimer geometry; therefore, a more effective linear hydrogen bond is obtained by pyramidalization of the amide group away from the aromatic rings, toward the oxygen of the second molecule. Such pyramidalization is observed in the benzamide43 (BZAMID02) and nicotinamide44 (NICOAM03) crystal structures, where neutron diffraction data were used to accurately determine the hydrogen atom positions. Unlike the results for forms III and IV, the NH2 geometry optimizations for forms I and II do not favor linear N-H- -O bonds. Instead, a lower total (conformational + lattice) energy


Cruz Cabeza et al.

Table 2. Structural Information and Calculated Energies of the Experimentally Observed and Corresponding Predicted Structures, Using the W99 + Atomic Multipoles Model Potential cryst struct

F (g/cm3)

a (Å)

b (Å)

c (Å)

exptl CBMZPN01 pred Hina pred Houtb pred min-confc planar, HNCN ) 0°

1.35 1.31 1.31 1.32

7.53 7.88 7.80 7.86

11.15 11.00 10.98 10.96

15.47 15.74 15.78 15.73

Form III, Space Group P21/c 90.0 116.17 90.0 90.0 118.73 90.0 90.0 117.97 90.0 90.0 118.45 90.0

exptl CBMZPN12 pred Hina pred Houtb pred min-confc Hout HNCN ) 14.6°

1.30 1.24 1.24 1.24

26.61 26.87 25.73 25.73

6.93 7.00 6.77 6.77

13.96 14.46 14.58 14.58

Form IV, Space Group C2/c 90.0 109.70 90.0 90.0 111.87 90.0 90.0 100.03 90.0 90.0 100.03 90.0

exptl CBMZPN03 pred Hina pred Houtb pred min-confc HNCN ) 30°

1.24 1.21 1.22 1.23

35.45 36.84 36.54 36.51

35.45 36.84 36.54 36.51

5.25 4.98 5.00 4.98

Form II, Space Group R3h 90.0 90.0 120.0 90.0 90.0 120.0 90.0 90.0 120.0 90.0 90.0 120.0

exptl CBMZPN11 pred Hina pred Houtb pred min-conf c HNCN ) 30°

1.34 1.27 1.27 1.29

20.57 21.46 21.41 21.33

5.17 5.01 5.01 5.00

22.24 23.13 23.20 23.03

R (deg)

β (deg)

γ (deg)

Form I, Space Group P1h 88.01 84.12 85.19 90.42 84.93 86.71 90.31 85.23 85.27 87.37 84.91 86.70

T (K)

∆E(conf) (kJ/mol)

E(lattice) (kJ/mol)

E(lattice) + ∆E(conf) (kJ/mol)

0 +1.96 +1.86

-113.62 -116.19 -116.77

-113.62 -114.23 -114.91

0 +1.96 +1.96

-112.84 -116.68 -116.68

-112.84 -114.72 -114.72

0 +1.96 +3.18

-104.74 -107.53 -109.71

-104.74 -105.57 -106.53

0 +1.96 +3.18

-108.14 -110.22 -111.87

-108.14 -108.26 -108.69

298

158

298

158

a Lattice parameters of the crystal structure predicted with the H molecular model. b Lattice parameters of the crystal structure predicted with the H in out molecular model. c Lattice parameters of the final predicted crystal structure, after determining the preferred amide geometry.

Figure 9. Geometries of the dimer motifs found in carbamazepine forms II-IV experimental (X-ray) crystal structures, distances between the NCO planes (d), and predicted conformations.

is obtained with amide pyramidalization away from the carbonyl oxygen of the second molecule; a negative HNCN angle is needed for a linear hydrogen bond, in contrast to that predicted from our calculations. Analysis of the energy-minimized structure suggests that this counterintuitive NH2 distortion is a result of stabilizing CH- - -N contacts with the aromatic ring of a third molecule. Indeed, such contacts (dH--N ) 2.57 Å) are increasingly considered to be weak hydrogen bonds,45 with structure-directing capabilities.

Analysis of the Hydrogen Bonding and Packing of the Structures. Although all experimental polymorphs exhibit the same hydrogen-bonding arrangement of anti-carboxamide dimers, other stable patterns are present in the energetically competitive predicted structures. The low-energy hypothetical crystal structures were analyzed and grouped by type of hydrogen-bonding motif (dimers, chains, and structures lacking typical hydrogen bonds): see Figure 8. Non-hydrogen-bonded hypothetical crystals are the least stable, though they represent some of the



Figure 10. Common aromatic patterns found in the observed crystals: (top) offset π stacking packing in the trigonal and triclinic polymorphs; (bottom) sandwich-herringbone packing in both monoclinic crystals as well as in the most stable but as yet unobserved chain structure.

densest structuressclearly, optimization of hydrogen-bonding interactions results in less efficient packing. It is interesting to note that hydrogen-bonded chains are as abundant and stable as dimers, although no structures with chain motifs have been experimentally observed. Only one of the three most stable structures exhibits chains, but there are many others more stable than both forms I and II. The lack of observed chain structures, even after an extensive polymorph screen,14 suggests the importance of molecular recognition effects during crystallization. Dimers pack as a unit and, if they are present in the solution or melt, chain formation would require disruption of this supramolecular synthon during crystal nucleation and growth. While such kinetic factors may be responsible for the lack of observed hydrogen-bonded chains, we note that the structurally similar 10,11-dihydrocarbamazepine46 (CSD refcode VACTAU01) does crystallize with the amide groups forming H-bonded chains. That such a small perturbation of the molecular structure switches the hydrogen bonding from dimer to chain shows that the dimer motif may not be so robust and further investigation of crystal growth factors, as well as entropic contributions to the crystal energies, is needed. Examination of the four observed structures indicates that it is the packing of the dimer unit47 that differentiates between

the known polymorphs. It is well-known that aromatic rings in crystals tend to organize themselves in one of two ways: in either a stacked or a herringbone fashion. In the stacking arrangement, the aromatic rings are aligned parallel (dz ) 3.13.4 Å) and shifted relative to each other (dxy ) 3.6-4.3 Å). In the herringbone arrangement, the preferred geometry has angles between the rings in the range of 40-90°, with the closest contacts between carbon and hydrogen atoms.48 We analyzed the computer-generated crystal structures for common patterns of aromatic interactions. The two most abundant are π-π stacking and edge-to-face interactions that were classified as offset π stacking and sandwich-herringbone (Figure 10);49,50 the three most stable crystals showed sandwich-herringbone packing, while both forms I and II, as well as most of the other low-energy predictions, showed offset π stacking. Crystal structures are a balance of many types of interactions, and while hydrogen bonding is very important, the aromatic interactions also appear to be important in this case and separate the most stable from the metastable crystals. The majority of the carbamazepine molecular surface is nonpolar, and interactions between these hydrophobic surfaces must contribute importantly to the packing energy, with perhaps as much structure-directing influence as the hydrogen bonds. We are


exploring the importance of the non-hydrogen-bonding interactions in carbamazepine and related molecules in a separate study. Conclusions The pharmaceutical compound carbamazepine exhibits a certain degree of conformational flexibility, due to the pyramidalization of the amide group nitrogen. We have demonstrated that the influence of this apparently small amount of flexibility on lattice energy calculations plays a significant role in the modeling and crystal structure prediction of the polymorphs of carbamazepine. Many of the low-energy predicted crystal structures gave different optimal amide geometries, and the relaxation of this molecular degree of freedom has an important effect on the energy ranking of the predicted structures. The calculated relative stabilities of the polymorphs show qualitative agreement with the available experimental datas only the stability of form IV was in disagreement with the one experimental study available for this polymorph. The final calculations found the two known monoclinic polymorphs (forms III and IV) as the most stable possible structures, together with a third, as yet unobserved, crystal that contains the alternative chain hydrogen-bonding motif. The remaining polymorphs (forms I and II) were found at 6 and 8 kJ/mol above the global minimum, higher in energy than many computergenerated crystal structures; these polymorphs would not have been suggested as likely observed structures by lattice energy calculations alone. Entropic or crystal growth considerations are suggested therefore as important factors to be considered in any further analysis. Acknowledgment. We thank Dr. Neil Feeder and Dr. Pete Marshall for valuable discussions regarding this project. We also thank Prof. Matzger for providing us with the structure of carbamazepine form IV, Dr. James Chisholm for the COMPACK algorithm implemented in COSET, and one of the reviewers for drawing our attention to ref 47. Finally, we thank the Pfizer Institute for Pharmaceutical Materials Science for funding. Supporting Information Available: Text and a figure detailing the electrostatic model dependence of the conformation optimizations and a CIF file giving the 10 lowest energy structures from the crystal structure prediction calculations. This material is available free of charge via the Internet at http://pubs.acs.org.

References (1) Bernstein, J. Polymorphism in Molecular Crystals; Oxford University Press: Oxford, U.K., 2002. (2) Bauer, J.; Spanton, S.; Henry, R.; Quick, J.; Dziki, W.; Porter, W.; Morris, J. Pharm. Res. 2001, 18, 859-866. (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) Price, S. L. AdV. Drug DeliV. ReV. 2004, 56, 301-319. (5) Day, G. M.; Motherwell, S.; Ammon, H. L.; Boerrigter, S. X. M.; Della Valle, R. G.; Venuti, E.; Dzyabchenko, A.; Dunitz, J. D.; Schweizer, B.; van Eijck, B. P.; Facelli, J. C.; Bazterra, V. E.; Ferraro, M. B.; Hofmann, D. W. M.; Leusen, F. J. J.; Liang, C.; Pantelides, C. C.; Karamertzanis, P. G.; Price, S. L.; Lewis, T. C.; Nowell, H.; Torrisi, A.; Scheraga, H.; Arnautova, Y. A.; Schmidt, M. U.; Verwer, P. Acta Crystallogr., Sect. B: Struct. Sci. 2005, 61, 511-527. (6) Mooij, W. T. M.; Leusen, F. J. J. Phys. Chem. Chem. Phys. 2001, 3, 5063-5066. (7) 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.;

Cruz Cabeza et al.

(8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) (25) (26)

(27) (28) (29) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39) (40) (41)

(42) (43) (44) (45) (46) (47) (48) (49) (50)

Williams, D. E. Acta Crystallogr., Sect. B: Struct. Sci. 2002, 58, 647-661. Ouvrard, C.; Price, S. L. Cryst. Growth Des. 2004, 4, 1119-1127. Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Am. Chem. Soc. 2000, 122, 3500-3505. Reboul, J. P.; Cristau, B.; Soyfer, J. C.; Astier, J. P. Acta Crystallogr., Sect. B: Struct. Sci. 1981, 37, 1844-1848. Himes, V. L.; Mighell, A. D.; De Camp, W. H. Acta Crystallogr., Sect. B: Struct. Sci. 1981, 37, 2242-2245. Lowes, M. M. J.; Caira, M. R.; Lo¨tter, A. P.; Van Der Watt, J. G. J. Pharm. Sci. 1987, 76, 744-752. Grzesiak, A. L.; Lang, M. D.; Kim, K.; Matzger, A. J. J. Pharm. Sci. 2003, 92, 2260-2271. Florence, A. J.; Johnston, A.; Price, S. L.; Nowell, H.; Kennedy, A. R.; Shankland, N. J. Pharm. Sci., in press. Mannfors, B. E.; Mirkin, N. G.; Palmo, K.; Krimm, S. J. Phys. Chem. A 2003, 107, 1825-1832. Mannfors, B.; Palmo, K.; Krimm, S. J. Mol. Struct. 2000, 556, 1-21. Momany, F. A. J. Phys. Chem. 1978, 82, 592. Cox, S. R.; Williams, D. E. J. Comput. Chem. 1981, 2, 304. Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047-1064. MS Modelling, Release 3.0.1; Accelrys Inc., San Diego, CA, 2004. Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13288-13249. Delley, B. J. J. Chem. Phys. 1990, 92, 508-517. Williams, D. E. J. Mol. Struct. 1999, 486, 321-347. Williams, D. E. J. Comput. Chem. 2001, 22, 1154-1166. Williams, D. E. J. Comput. Chem. 2001, 22, 1-20. Amos, R. D. CADPAC, version 6.0; 1995 (with contributions from Alberts, I. L.; Andrews, J. S.; Colwell, S. M.; Handy, N. C.; Jayatilaka, D.; Knowles, P. J.; Kobayashi, R.; Koga, N.; Laidig, K. E.; Maslen, P. E.; Murray, C. W.; Rice, J. E.; Sanz, J.; Simandiras, E. D.; Stone, A. J.; Su, M. D.). Karfunkel, H. R.; Gdanitz, R. J. J. Comput. Chem. 1992, 13, 1171. Karfunkel, H. R.; Leusen, F. J. J. Speedup 1992, 6, 43. Verwer, P.; Leusen, F. J. J. ReV. Comput. Chem. 1998, 12, 327365. Cerius2, version 4.6; Accelrys Inc., San Diego, CA, 1997. Day, G. M.; Chisholm, J.; Shan, N.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2004, 4, 1327-1340. Price, S. L.; Willock, D. J.; Leslie, M.; Day, G. M. DMAREL, version 3.1; 2001. Chisholm, J. A.; Motherwell, S. J. Appl. Crystallogr. 2005, 38, 228231. Day, G. M.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2005, 5, 1023-1033. Edwards, A. D.; Shekunov, B. Y.; Forbes, R. T.; Grossmann, J. G.; York, P. J. Pharm. Sci. 2001, 90, 1106-1114. Behme, R. J.; Brooke, D. J. Pharm. Sci. 1991, 80, 986-990. Urakami, K.; Shono, Y.; Higashi, A.; Umemoto, K.; Godo, M. Chem. Pharm. Bull. 2002, 50, 263-267. Park, K.; Evans, J. M. B.; Myerson, A. S. Cryst. Growth Des. 2003, 3, 991-995. Gu, C. H.; Grant, D. J. W. J. Pharm. Sci. 2001, 90, 1277-1287. Ceolin, R. T. S.; Gardette, M. F.; Agafonov, V. N.; Dzyabchenko, A. V.; Bachet, B. J. Pharm. Sci. 1997, 86, 1062-1065. Harris, R. K.; Ghi, P. Y.; Puschmann, H.; Apperley, D. C.; Griesser, U. J.; Hammond, R. B.; Ma, C.; Roberts, K. J.; Pearce, G. J.; Yates, J. R.; Pickard, C. J. Org. Process Res. DeV. 2005, 9, 902-910. Gavezzotti, A.; Filippini, G. J. Am. Chem. Soc. 1995, 117, 1229912305. Gao, Q.; Jeffrey, G. A.; Ruble, J. R.; McMullan, R. K. Acta Crystallogr., Sect. B: Struct. Sci. 1991, 47, 742. Miwa, Y.; Mizuno, T.; Tsuchida, K.; Taga, T.; Iwata, Y. Acta Crystallogr., Sect. B: Struct. Sci. 1999, 55, 78. Allen, F. H.; Bird, C. M.; Rowland, R. S. Acta Crystallogr., Sect. B: Struct. Sci. 1995, B51, 1068-1081. Bandoli, G.; Nicolini, M.; Ongaro, A.; Volpe, G.; Rubello, A. J. Cryst. Spectrosc. Res. 1992, 177. Gelbrich, T.; Hursthouse, M. B. CrystEngComm 2006, 8, 449-461. Hunter, C. A.; Sanders, J. K. M. J. Am. Chem. Soc. 1990, 112, 55255534. Gavezzotti, A. Acta Crystallogr., Sect. B: Struct. Sci. 1988, 44, 427434. Desiraju, G. R.; Gavezzotti, A. Acta Crystallogr., Sect. B: Struct. Sci. 1989, 45, 473-482.

CG0601756

Amide Pyramidalization in Carbamazepine - American Chemical Society

Recommend Documents