Predicting the Spontaneous Chiral Resolution by Crystallization of a

Jun 28, 2008 - Matthew D. Gourlay, John Kendrick, and Frank J. J. Leusen*. Institute of Pharmaceutical InnoVation, UniVersity of Bradford, Bradford BD...
0 downloads 0 Views 802KB Size
Predicting the Spontaneous Chiral Resolution by Crystallization of a Pair of Flexible Nitroxide Radicals Matthew D. Gourlay, John Kendrick, and Frank J. J. Leusen* Institute of Pharmaceutical InnoVation, UniVersity of Bradford, Bradford BD7 1DP, United Kingdom

CRYSTAL GROWTH & DESIGN 2008 VOL. 8, NO. 8 2899–2905

ReceiVed December 20, 2007; ReVised Manuscript ReceiVed April 8, 2008

ABSTRACT: The separation of racemates into pure enantiomers through crystallization is an important industrial process. This study provides further validation of a novel, predictive approach for spontaneous resolution in which crystal structure prediction simulations are used to explore the relative stabilities of racemic solids versus enantiopure solids. 2-(4-Hydroxyphenyl)-2,5,5trimethylpyrrolidine-1-oxy (compound 1) has previously been shown to be a racemic conglomerate, while a similar compound, 2-(3-hydroxyphenyl)-2,5,5-trimethylpyrrolidine-1-oxy (compound 2), was not. A conformational search using the Dreiding force field revealed 10 conformational minima for compound 1, and 20 for compound 2. Atomic charges were calculated using unrestricted DFT B3LYP 6-311G** optimized structures, and a crystal structure prediction was performed using the Dreiding force field, considering all low-energy gas-phase conformations and all relevant space groups. Analysis of the predicted crystal structures suggests that compound 1 is a racemic conglomerate, but compound 2 is not. This is in agreement with the experimental evidence. Introduction The makeup of the human body is innately chiral, a fact which often requires the pharmaceutical industry to manufacture pure enantiomers of chiral drugs. While there are several routes to their manufacture, one of the most cost-effective is resolution by crystallization.1 The ability to predict whether a racemic mixture of a chiral molecule will crystallize as enantiopure crystals, or not, would therefore be of great economic benefit. In a previous paper, crystal structure prediction methods were used to study the relative stabilities of enantiopure and racemic crystals of 4-hydroxymethyl-2-oxazolidinone and 5-hydroxymethyl-2-oxazolidinone.2 The results indicated that the approach is promising, but further validation work was still required. In the present work, this application is extended to two related flexible molecules containing a nitroxide radical group; 2-(4hydroxyphenyl)-2,5,5-trimethylpyrrolidine-1-oxy (compound 1) and 2-(3-hydroxyphenyl)-2,5,5-trimethylpyrrolidine-1-oxy (compound 2), shown in Figure 1 along with a definition of the important geometric variables. As in the previous work, the two molecules are structurally similar (differing only in the position of the hydroxyl group on the phenol ring), but they exhibit different crystallization behavior. These compounds were first crystallized to investigate the use of paramagnetic materials in liquid crystals.3 Using IR, ESR, and crystallization, it has been shown that compound 1 forms a racemic conglomerate, while compound 2 does not.3,4 The radical N-O group in these compounds is expected to present some challenges, as the commonly used molecular mechanics force fields have not been parametrized for radicals. In the past, radical systems have mostly been investigated using open shell ab initio quantum mechanics (QM) techniques.5,6 There has been a previous attempt at crystal structure prediction of a radical system using PROMET, where a rigid molecule, 2-hydronitronylnitroxide, was investigated.7 In this study, crystal structure prediction methods are used to predict whether racemic crystals or enantiopure crystals will precipitate from a racemic solution of compound 1 or compound 2. The results are then compared with the experimental evidence. * Author to whom correspondence should be addressed. Tel: +44(0)1274 236144. Fax: +44(0)1274 236155. E-mail: [email protected].

The use of lattice energies to predict the outcome of crystallization8 ignores the effects of kinetics and of solvent. This assumption is supported by the conclusions of a previous study looking at diastereomeric salt resolutions.8,9 Methodology Molecular mechanics calculations were performed using Cerius2 v 4.610 and Materials Studio v 4.0.11 Ab initio calculations were performed using GAMESS-UK v 6.3.12 The Cambridge Structural Database (5.28, November 2006)13 (CSD) was searched to find the crystal structures of the compounds. For all solid-state simulations, the Ewald method14,15 was used for the electrostatic and van der Waals interaction terms. Gasteiger charges were used for an initial conformational search. Subsequent calculations used electrostatic potential derived charges. As the crystal structure prediction method uses a rigid body approximation in the initial search for crystal packing alternatives, it is necessary to perform a conformational analysis to determine all lowenergy conformations to be used as input for the packing calculations. The molecules were drawn and conformational grid scans were performed using Cerius2.10 The force field used was Dreiding 2.2116 with Gasteiger charges (as implemented in the Cerius2 and Materials Studio packages). Initially a short molecular dynamics simulation was performed using the Forcite module of Materials Studio (simulation length of 6 ps with a 1 fs time step at a temperature of 600 K, taking conformations every 100 steps for an independent optimization) to assess the flexibility of the ring. A grid scan was performed on all conformations with a distinct ring geometry by varying torsion angles T1 and T2 between 0° and 360° in 10° increments. The optimized lowenergy conformations were reoptimized by GAMESS-UK12 using an unrestricted B3LYP density functional17,18 with the 6-311G** basis set.19,20 In contrast to the previous study,2 where the point charges were calculated from the QM electron distributions for each conformation independently and then averaged over the conformations using Boltzmann weighting, in this study the electrostatic potential was averaged over all conformations using Boltzmann weighting before calculating the point charges from the averaged electrostatic potential. Using these atomic charges, all low-energy conformations were reoptimized using the Dreiding 2.21 force field. These optimized gasphase conformations were used as the starting points for crystal structure prediction using the Materials Studio polymorph predictor.11 The polymorph predictor was set to its default fine setting; however, the packing step was altered to its medium setting (this sets the simulated annealing algorithm to a temperature range of 300-60000 K with a cooling factor of 0.002, requiring 10 consecutive steps to be accepted before cooling and a maximum of 4000 steps). The 13 most

10.1021/cg701256e CCC: $40.75  2008 American Chemical Society Published on Web 06/28/2008

2900 Crystal Growth & Design, Vol. 8, No. 8, 2008

Gourlay et al.

Figure 1. The following notation is used to describe the torsion angles in the molecules. Ring puckering is described by torsion angles RT1 (N5-C4-C3-C2) and RT2 (C4-C3-C2-C1). The pyramidal nature of the NO group is defined by an improper torsion ITNO (C4-C1-N5-O6); here an angle of 180° means the group is flat. The rotation of the phenol ring is defined by torsion angle T1 (N5-C4-C10-C15). The rotation of the hydroxyl group is defined by torsion angle T2 (C14-C13-016-H17 in compound 1 and C15-C14-O16-H17 in compound 2). An asterisk indicates a chiral center. Table 1. Analysis of the Effects of Different Atom Types on the Geometry of the N-O Moietya structure

atom typesb

EYABAH (compound 1) WASHII (compound 2) compound 1 N_2, O_3 compound 1 N_2, O_2 compound 1 N_3, O_3 compound 2 N_2, O_3 compound 2 N_2, O_2 compound 2 N_3, O_3

N-O bond representation

double double single double double single

N-O length (Å)

ITNO (°)

1.28 1.28 1.28 1.17 1.36 1.28 1.17 1.36

170.6 179.2 178.8 178.8 145.7 178.5 179.1 155.6

a The first two rows show the experimental values. The remainder shows results for molecules minimized in the gas-phase using the Dreiding force field and Gasteiger charges. The rows highlighted in bold indicate the atom types chosen for the force field representation used in the rest of the paper. b In the Dreiding force field _2 and _3 refer to a trigonal (sp2) and a tetrahedral (sp3) atom, respectively.

common space groups21 found in organic crystals registered in the CSD were selected (P21/c, P1j, P212121, P21, C2/c, Pbca, Pnma, Pna21, Pbcn, P1, Cc, C2, and Pca21). Clustering of the predicted polymorphs was done using the polymorph clustering routine in Materials Studio, combined with a spreadsheet macro. The Materials Studio clustering routine calculates a crystal similarity measure (CSM) from the list of interatomic distances for each combination of atoms or force field atom types in the crystal. The CSM is a number between 0 and 1. A value of 0 indicates that the two crystals are identical, while 1 implies no crystal similarity. The Materials Studio polymorph clustering routine was unable to cluster all the polymorphs generated in this investigation within a reasonable amount of disk space and time. Therefore, the routine was only used to cluster all the polymorphs for each conformation used in the prediction. An Excel macro was written which compares the lattice energy, density, and the CSM for each crystal structure, and this macro was subsequently used to cluster the unique polymorphs obtained for the different conformations. The polymorph predictions were repeated five times for each starting conformation to ensure adequate sampling of the crystal configurations.

Results and Discussion A search in the CSD13 found one crystal structure for compound 1 (reference code EYABAH) and one crystal structure for compound 2 (reference code WASHII). Table 1 shows a summary of the geometry of the NO group obtained by minimizing the energy of the molecules in the gas phase, using different Dreiding atom types. The experimental N-O · bond lengths and ITNO angles were used to check the suitability of the force field. Unfortunately, there is no direct way to represent the N-O · radical in the Dreiding force field

Figure 2. Charge assignments for the phenyl ring of compound 1 based on the ESP of different conformations. The charges in panels 2A and 2B are determined by fitting to the ESP of the conformations shown. Panel 2C shows the results of fitting the charges to a Boltzmann weighted average of the ESPs for all conformations.

as it has not been parametrized for this class of compounds. However, the N-O · radical is expected to be stabilized by resonance with N+•-O-.22,23 An acceptable representation of the radical structures was found to be N+dO-, where the + and - indicate the formal charges placed on the atoms. Using this representation, indicated by Dreiding force field atom types N_2 and O_3 in Table 1, the bond length and ITNO angle are close to the experimental values. No further changes were made to the force field. Two distinctive ring conformations were found in the molecular dynamics simulation, one with a positive RT1 value, the other with a negative RT1 value. The two lowest energy conformers, representative of the different ring conformations, found during the molecular dynamics simulation were geometry optimized and used as the starting point for a detailed conformational search for each compound using a grid scan method with Gasteiger charges. This resulted in 10 unique, minimized conformations for compound 1 and 20 for compound 2. The contour plots of energy against torsion angles T1 and T2 (shown in the Supporting Information) show four energy wells for each

Predicting Spontaneous Chiral Resolution

Crystal Growth & Design, Vol. 8, No. 8, 2008 2901

Table 2. Deviations of the Calculated Lattice Parameters from Experiment, Using Different Charge Setsa unit cell dimensions charge set

N-O bond length (Å)

b (Å)

c (Å)

β (°)

% RMS (abc)

EYABAH compound 1 compound 1

Gasteiger B3LYP 6-311G**

1.28 1.28 1.30

8.614 0.04% 0.92%

10.721 1.49% 0.58%

13.421 0.53% 0.72%

90.0 0.0% 0.0%

0.00 0.91 0.75

WASHII compound 2 compound 2

Gasteiger B3LYP 6-311G**

1.28 1.28 1.30

8.122 0.42% -1.48%

20.216 5.76% 4.87%

8.094 3.66% -0.08%

110.4 0.5% -1.3%

0.00 3.95 2.94

a (Å)

a The calculations used the Dreiding force field and the N+dO- representation of the N-O bond. % RMS (abc) is the root-mean-squared value of the percentage deviation values for the lattice lengths.

ring conformation. Each well is flanked by shoulders that are slightly higher in energy. Note that the number of unique conformations found for compound 1 is half the number found for compound 2. This is because of the torsional equivalence found in the phenol substituent of compound 1; namely, simultaneous rotation by 180° around T1 and T2 does not change the molecular conformation. The conformations found by the search were optimized using UDFT (B3LYP, 6-311G**) calculations. The number of unique conformations was reduced to 4 for compound 1 and 8 for compound 2. Considering the conformational landscapes, it is not surprising that the DFT calculations found fewer conformations. The higher energy shoulders are probably an artifact of the molecular mechanics force field. This was not confirmed with the QM methods as a full conformational search with QM was not feasible. Figure 2, panels A and B, shows that the charge assignments based on the electrostatic potential (ESP) for atoms C12 and C14 of compound 1 depend strongly on the orientation of the hydroxyl group. Ideally, this should not be the case as C12 and C14 are torsionally equivalent. For torsional equivalence the molecule maintains its geometry and electrostatic distribution after rotation around a given torsion angle. In this work equivalence has not been imposed, as the fit of the resulting charges to the ESP is poorer. Instead the ESPs of all the unique, QM optimized conformers were included in the fitting of the atomic charges.24 A Boltzmann factor at 298 K was used to weight the ESP of each conformer according to its energy relative to the global minimum. The final charge assigment for compound 1 is shown in Figure 2, panel C, while the complete charge assignment for both compounds is available in the Supporting Information. Although carbons C11 and C15 (and the hydrogens attached to these carbons) in compound 1 are also torsionally equivalent, assignment of the ESP charges means that their electrostatic properties are not equivalent due to their different interaction with the nitroxyl group. Simultaneous 180° rotation about T1 and T2 followed by optimization should not alter the molecular geometry in the crystal, but due to these electrostatic differences, the force field calculation produces slightly different lattice energies and geometries. This problem is addressed by discarding any predicted crystal structure where T1 is outside the range -90 to +90°. This is not an issue for compound 2, as its phenol group has a meta substituted hydroxyl group, and none of its atoms are torsionally equivalent. To test the suitability of the molecular mechanics model, the experimentally determined crystal structures were geometry optimized. Table 2 shows that the Gasteiger charges which were used for the initial conformational analysis result in larger deviations between the experimental and force field optimized lattice parameters than the potential derived charges. Subsequent calculations were therefore performed with the ESP charges.

Table 3. Torsion Angles and Relative Energies of the Lowest Energy Conformers of Compound 1 Obtained with the Dreiding Force Field and ESP Chargesa conformation

ITNO (°)

RT1 (°)

RT2 (°)

T1 (°)

T2 (°)

∆ energy (kcal mol-1)

Crystal_SM Conf 1 Conf 2 Conf 3 Conf 4 Conf 5 Conf 6 Conf 7 Conf 8 Conf 9 Conf 10

179.4 178.7 178.7 179.4 179.4 -178.3 -178.3 177.2 177.2 176.8 176.8

-28.0 19.6 19.5 -28.0 -28.0 -36.6 -36.6 -13.7 -13.7 27.1 27.1

32.6 -28.9 -28.9 32.6 32.6 33.7 33.7 23.2 23.3 -32.2 -32.2

-4.7 0.0 0.0 -4.7 -4.7 -58.5 -58.4 49.5 48.9 -50.9 -50.7

-0.1 -0.3 -180.0 -0.1 -179.9 0.0 -180.0 -179.8 -0.2 0.4 -179.8

0.07 0.00 0.05 0.07 0.11 2.10 2.19 3.08 3.11 3.69 3.79

a Crystal_SM represents a single molecule extracted from the experimental crystal structure that has been optimized. Torsion angles are defined in Figure 1. The N-O bond length for all conformers was found to be constant at 1.29 Å. The line divides the four low-energy conformations from the higher energy conformations.

The ESP charges were applied to the minimum energy conformations generated by the grid scan. These conformations were reoptimized using the Dreiding force field and used as starting points for the crystal structure prediction simulations. The conformations are detailed in Tables 3 and 4. Analysis of the energies of the gas phase conformations reveals clusters of conformations with similar energies. For compound 1 (see Table 3), the low energy cluster contains four conformations within 0.11 kcal mol-1 of each other, followed by a 2 kcal mol-1 gap to the next cluster of higher energy conformations. For compound 2 (see Table 4), this clustering in conformational energies is less pronounced. Here, the first eight conformations are found within an energy window of 0.70 kcal mol-1, followed by a gap of 1.4 kcal mol-1 to the cluster of higher energy conformations. As the molecules are treated as completely flexible in the final lattice energy minimization step during the polymorph prediction approach used in this study, any high energy conformations may optimize to low energy conformations. This allows us to reduce the computational burden of the crystal structure prediction simulations. These simulations are performed several times to ensure a thorough sampling of the configurational space. In this study, the simulations were performed five times for the low energy conformations (conformations 1-4 of compound 1, and conformations 1-8 of compound 2), but only twice for the higher energy conformations. The latter runs were carried out to check whether any higher energy conformations may lead to stable unique crystal packing alternatives. Tables 5 and 6 show details of the predicted lowest energy enantiopure and racemic crystal structures for compounds 1 and 2, respectively. In the 100 lowest energy structures for compound 1, only one structure has been discarded due to rotation of T1 by more than 90°. The experimental crystal structure of

2902 Crystal Growth & Design, Vol. 8, No. 8, 2008

Gourlay et al.

Table 4. Torsion Angles and Relative Energies of the Lowest Energy Conformers of Compound 2 Obtained with the Dreiding Force Field and ESP Chargesa conformation

ITNO (°)

RT1 (°)

Crystal_SM Conf 1 Conf 2 Conf 3 Conf 4 Conf 5 Conf 6 Conf 7 Conf 8 Conf 9 Conf 10 Conf 11 Conf 12 Conf 13 Conf 14 Conf 15 Conf 16 Conf 17 Conf 18 Conf 19 Conf 20

179.4 179.4 178.7 179.2 179.2 178.7 179.2 178.6 178.6 -178.5 -178.6 -178.6 -178.6 177.0 177.1 177.1 177.2 176.7 177.0 176.9 176.7

-28.0 -28.0 19.6 -28.2 -27.9 19.3 -27.9 19.6 19.5 -36.3 -35.8 -36.0 -36.3 -13.7 -13.4 -13.6 -13.7 27.0 27.8 27.9 27.0

RT2 (°)

T1 (°)

∆ energy (kcal mol-1)

T2 (°)

32.7 -4.7 -0.6 32.7 -4.7 -0.6 -29.0 0.0 -0.8 32.7 -4.3 -179.7 32.7 175.9 179.9 -28.9 0.5 -179.8 32.7 176.0 0.2 -28.9 179.2 180.0 -28.9 179.3 0.3 33.8 -58.6 0.4 33.5 119.7 180.0 33.5 119.9 0.1 33.7 -58.8 179.9 23.0 51.3 -179.8 22.9 -129.7 0.4 23.1 -130.5 179.9 23.2 50.1 -0.5 -32.2 -51.0 0.7 -32.3 129.1 179.8 -32.3 129.4 -0.3 -32.3 -51.1 180.0

0.00 0.00 0.09 0.35 0.53 0.54 0.56 0.67 0.70 2.08 2.21 2.36 2.60 2.77 2.78 2.82 2.85 3.55 4.02 4.17 4.40

a Crystal_SM represents a single molecule extracted from the experimental crystal structure that has been optimized. Torsion angles are defined in Figure 1. The N-O bond length for all conformers was found to be constant at 1.29 Å. The line divides the eight low-energy conformations from the higher energy conformations.

compound 1 (enantiopure structure 2 in Table 5; see also Figure 3) is 0.13 kcal mol-1 above the lowest energy structure. For compound 2, the experimental crystal structure is 0.88 kcal mol-1 above the lowest energy racemic structure (racemic structure 10 in Table 6; see also Figure 4). The calculations correctly predict an enantiopure crystal structure to have the lowest lattice energy for compound 1, and a racemic crystal structure for compound 2. The predictions show an energy gap of 0.57 kcal mol-1 between the lowest energy racemic and enantiopure crystal structures generated for compound 2. For

compound 1, the racemic conglomerate, this energy gap is smaller (0.19 kcal mol-1). Tables 7 and 8 show more information on the molecular geometries in the predicted crystal structures. The N-O bond length is found to be constant at 1.30 Å, and the ITNO torsion angle, in most cases, remains around 180°. It is noteworthy that the geometries of both molecules in the various predicted crystal structures do not deviate significantly from the gas phase geometries found in the conformational analysis. Intermolecular hydrogen bonding interactions between the hydroxyl and nitroxide groups are observed in all the crystal structures. Despite the potential for OH-OH hydrogen bonding, no such bonding is seen. The formation of hydrogen bonding chains is the most common pattern (an example is shown in Figure 3). Hydrogen bonding chains propagate infinitely in one direction through the crystal. RS dimers are found only in racemic crystal structures and are observed experimentally for compound 2, as can be seen in Figure 4. This implies that for compound 2, two molecules of opposite chirality pack more efficiently than two molecules of the same chirality. Considering all the crystal structures generated, there is a higher percentage of RS dimer structures for compound 2 than for compound 1. The force field dictates both the configurations of the molecules found in the crystal and their relative energies. The former reflects the presence of a local potential energy well, and even relatively poor force fields can predict the presence of such minima, without necessarily accurately predicting their relative energies. Thus, the observation of a common structural feature such as the RS dimer may be an indicator that the system will form a racemic conglomerate. The assumption that all relevant crystal structures can be generated starting from a subset of the complete set of conformations holds true, as shown by the last columns in Tables 7 and 8. Further analysis of the 30 lowest energy enantiopure and racemic structures for both compounds (of which only the first 10 are shown in Tables 7 and 8) shows that they can all be

Table 5. The 10 Lowest-Energy Predicted Enantiopure and Racemic Crystal Structures of Compound 1a compound 1 racemic

space group

1 2 3 4 5 6 7 8 9 10

Pbca P21/c Pbca P21/c P21/c P21/c Pna21 Pna21 P21/c Pbca

compound 1 enantiopure

space group

1 2 3 4 5 6 7 8 9 10

P212121 P212121 P21 P212121 P212121 P212121 P212121 P212121 P21 P21

density (g cm-3)

∆ energy (kcal mol-1)

a (Å)

b (Å)

c (Å)

R (°)

β (°)

γ (°)

CSM

1.157 1.208 1.150 1.145 1.111 1.177 1.145 1.136 1.129 1.115

0.19 0.65 0.77 0.82 0.99 1.01 1.01 1.11 1.21 1.23

11.80 15.08 16.66 9.98 9.86 9.81 13.29 13.54 7.47 14.15

13.53 11.44 11.20 10.87 14.63 13.59 10.49 9.62 11.46 14.61

15.83 14.36 13.63 14.37 16.28 12.74 9.16 9.89 16.32 12.70

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

90.0 150.7 90.0 125.0 145.9 133.0 90.0 90.0 112.0 90.0

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

0.195 0.247 0.218 0.230 0.213 0.216 0.206 0.193 0.243 0.216

density (g cm-3)

∆ energy (kcal mol-1)

a (Å)

b (Å)

c (Å)

R (°)

β (°)

γ (°)

CSM

1.161 1.155 1.142 1.154 1.130 1.179 1.140 1.111 1.123 1.096

0.00 0.13 0.76 1.09 1.21 1.56 1.98 2.06 2.10 2.22

13.84 13.52 8.40 13.34 10.08 6.96 7.37 7.23 7.53 9.89

8.50 10.78 10.27 10.75 12.54 17.54 11.19 15.38 11.52 11.63

10.71 8.69 7.60 8.84 10.24 10.16 15.55 11.84 8.13 6.88

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

90.0 90.0 77.5 90.0 90.0 90.0 90.0 90.0 112.6 122.5

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

0.216 0.018 0.223 0.191 0.231 0.194 0.231 0.199 0.205 0.214

a The relative lattice energy is reported in kcal mol-1 per asymmetric unit. The Crystal Similarity Measure (CSM) is used with the optimized experimental crystal structure as the reference structure; a value close to 0 indicates a match between the reference and tested crystals. In this case, enantiopure structure 2 matches the experimental crystal structure and is highlighted in bold.

Predicting Spontaneous Chiral Resolution

Crystal Growth & Design, Vol. 8, No. 8, 2008 2903

Table 6. The 10 Lowest-Energy Predicted Enantiopure and Racemic Crystal Structures of Compound 2a compound 2 racemic

space group

1 2 3 4 5 6 7 8 9 10

P21/c Pbca P21/c P21/c C2/c P21/c P1j P21/c P21/c P21/c

compound 2 enantiopure

space group

1 2 3 4 5 6 7 8 9 10

P212121 P212121 P212121 P212121 P21 P212121 P212121 P212121 P21 P212121

density (g cm-3)

∆ energy (kcal mol-1)

a (Å)

b (Å)

c (Å)

R (°)

β (°)

γ (°)

CSM

1.170 1.149 1.156 1.140 1.174 1.156 1.160 1.125 1.149 1.128

0.00 0.14 0.29 0.49 0.62 0.67 0.74 0.79 0.80 0.88

9.29 11.26 7.98 13.20 12.21 17.94 7.77 12.42 14.42 8.00

11.46 13.58 11.32 8.45 9.52 13.87 9.47 8.44 10.76 21.20

13.82 16.65 14.26 18.14 21.46 13.37 13.39 16.96 13.82 8.09

90.0 90.0 90.0 90.0 90.0 90.0 117.1 90.0 90.0 90.0

58.2 90.0 100.8 140.7 88.5 157.7 92.7 47.0 143.6 71.0

90.0 90.0 90.0 90.0 90.0 90.0 126.7 90.0 90.0 90.0

0.250 0.217 0.208 0.226 0.209 0.249 0.225 0.210 0.235 0.006

density (g cm-3)

∆ energy (kcal mol-1)

a (Å)

b (Å)

c (Å)

R (°)

β (°)

γ (°)

CSM

1.139 1.121 1.138 1.148 1.129 1.150 1.100 1.131 1.137 1.099

0.57 1.27 1.46 1.55 1.56 1.61 1.77 1.82 1.83 1.84

13.98 8.95 12.54 12.73 7.60 12.72 14.28 8.47 8.19 13.29

8.59 16.49 9.98 9.94 11.39 6.93 11.04 13.97 11.57 10.73

10.69 8.85 10.27 10.08 8.05 14.44 8.44 10.94 7.12 9.33

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

90.0 90.0 90.0 90.0 68.5 90.0 90.0 90.0 107.3 90.0

90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0

0.214 0.223 0.264 0.226 0.251 0.241 0.212 0.235 0.225 0.234

a The relative lattice energy is reported in kcal mol-1 per asymmetric unit. The CSM is used with the optimized experimental crystal structure as the reference structure; a value close to 0 indicates a match between the reference and tested crystals. In this case, racemic structure 10 matches the experimental crystal structure and is highlighted in bold.

Figure 3. Predicted crystal structure for compound 1 (enantiopure structure 2 in Table 5), showing a hydrogen-bonded chain. This structure is available as a CIF file (Supporting Information).

generated starting from conformations 1-4 for compound 1 and conformations 1-8 for compound 2. Conclusions In this study, the ability of crystal structure prediction methods to calculate the relative stabilities of enantiopure and racemic crystals of 2-(4-hydroxyphenyl)-2,5,5-trimethylpyrrolidine-1-oxy (compound 1) and 2-(3-hydroxyphenyl)-2,5,5-trimethylpyrrolidine-1-oxy (compound 2) has been examined. The Dreiding force field was used with charges derived from the electrostatic potential calculated after unrestricted DFT (6-311G**) geometry optimizations. In the force field, the radical nature of these compounds was represented by setting the formal charges of

the oxygen atom to -1 and the nitrogen atom to +1, and selecting the corresponding force field atom types. This choice gives good predictions for the N-O bond length and the geometry around the nitrogen in the ring, as found in the CSD. The method used in this investigation predicts that the global minimum crystal structures belong to the same space groups as the experimentally determined crystal structures for both molecules. Unfortunately, the global minima do not correspond to the experimental crystals. For compound 1, the known crystal structure is predicted to be the second lowest energy crystal structure (0.1 kcal mol-1 above the global minimum) and for compound 2, the 11th predicted structure (0.9 kcal mol-1 above the global minimum) corresponds to the experimental crystal structure. Considering that the force field used for the lattice

2904 Crystal Growth & Design, Vol. 8, No. 8, 2008

Gourlay et al.

Figure 4. Predicted crystal structure for compound 2 (racemic structure 10 in Table 6), showing hydrogen-bonded RS dimers. This structure is available as a CIF file (Supporting Information). Table 7. The 10 Lowest-Energy Racemic and Enantiopure Crystal Structures for Compound 1 Showing the Conformational Details of the Asymmetric Unita compound 1 RT1 racemic (°) 1 2 3 4 5 6 7 8 9 10

-24.4 25.4 -34.8 -25.3 -25.5 -20.2 19.2 -20.4 22.7 24.3

compound 1 RT1 enantiopure (°) 1 2 3 4 5 6 7 8 9 10

19.1 -28.6 -23.9 25.9 20.8 27.4 27.2 -23.7 -26.9 24.4

RT2 (°)

T1 (°)

T2 (°)

31.9 2.7 -179.4 -29.6 18.4 -4.7 34.9 -14.3 -1.4 28.8 -3.4 -176.4 30.8 2.0 178.2 28.7 0.0 175.0 -25.7 2.6 -7.9 28.3 -2.5 -6.9 -30.7 3.4 -3.5 -26.6 14.9 -179.6

RT2 (°) -27.6 32.7 29.4 -31.0 -29.4 -30.4 -30.8 31.4 28.9 -29.8

T1 (°)

T2 (°)

1.5 4.7 -3.7 -0.2 4.4 -173.0 5.8 -2.4 1.7 -5.5 41.9 173.9 -0.2 179.6 0.9 -179.1 -3.6 -8.6 5.9 -171.1

hydrogen generated bonding from dimensionality conformations 1D 0D 1D 1D 1D 1D 1D 1D 1D 1D

4,6,7 1,9 3,5,8 4,6,7 4,6,7 4,6 1,9 3,5,8 1,9 2

hydrogen generated bonding from dimensionality conformations 1D 1D 1D 1D 1D 1D 1D 1D 1D 1D

1,8,9 3,5,8 4,6,7 1,9 1,9 2,10 2,10 4,6,7 3,5,8 2,10

Table 8. The 10 Lowest-Energy Racemic and Enantiopure Crystal Structures for Compound 2 Showing the Conformational Details of the Asymmetric Unita compound 2 RT1 racemic (°) 1 2 3 4 5 6 7 8 9 10

T1 (°)

T2 (°)

-25.0 31.5 -179.0 178.9 -34.2 34.8 -13.9 -178.4 -32.1 32.5 -8.7 -0.1 -28.3 32.2 -4.6 -0.7 -35.6 33.3 -13.8 -2.6 23.3 -28.0 -171.9 -3.1 -32.3 30.4 -10.6 -1.7 -25.6 30.4 -3.9 -2.4 -27.2 31.9 -4.3 -176.2 -32.9 33.0 -11.3 1.5

compound 2 RT1 enantiopure (°) 1 2 3 4 5 6 7 8 9 10

RT2 (°)

-28.5 -33.4 18.0 16.7 -28.5 23.2 20.6 18.2 -25.9 -25.7

RT2 (°)

T1 (°)

T2 (°)

32.7 -5.5 176.0 33.2 -13.0 1.8 -28.9 -3.4 -176.4 -28.3 -2.1 -174.3 31.4 176.0 -3.3 -28.2 5.1 -174.7 -28.1 4.1 175.8 -27.5 -1.0 177.1 32.3 178.0 172.2 30.8 -3.7 -174.7

hydrogen generated bonding from dimensionality conformations 1D 1D 0D 0D 0D 1D 0D 0D 1D 0D

4,10,15 3,12,13 1,9,16 1,11 1,9,16 8,19 1,2,9,16 1 3,12,13 1,9,16

hydrogen generated bonding from dimensionality conformations 1D 1D 1D 1D 1D 1D 1D 1D 1D 1D

3 1,9,16 5,20 5,20 6,11,14 5,20 5,20 5,20 4,10,15 3,12,13

a In the last but one column, 0D hydrogen bonding signifies dimer formation, while 1D signifies a chain. The final column shows the gas-phase conformations from which each structure was generated, as listed in Table 3.

a In the last but one column, 0D hydrogen bonding implies dimer formation, while 1D signifies a chain. The final column shows the gas-phase conformations from which each structure was generated, as listed in Table 4.

energy calculations in this study is not parametrized for these radical compounds, this energetic error of about 1 kcal mol-1 is acceptable. However, it is clear from the results that, in order to predict the correct experimental crystal structures and their relative energies, a method of calculating lattice energies with an error of much less than 1 kcal mol-1 is required. Recent progress in crystal structure prediction by Neumann et al.25 indicates that the use of solid-state density functional theory, coupled with an empirical correction for the van der Waals

interactions between atoms, is sufficient to provide the required accuracy. Such calculations are being performed on the crystal structures predicted in the present work and will be reported in due course. In this study it has been assumed that the low-energy molecular conformations found in the gas phase can be used as starting points for the generation of all low-energy crystal packing alternatives. Any distortions from the gas phase conformations due to intermolecular packing forces in the solid

Predicting Spontaneous Chiral Resolution

state are modeled by the force field. The assumption that all the predictable crystal structures can be generated from a selection of just the lower energy conformations has been shown to be valid, as consideration of the higher energy conformations produced no further stable crystal packing alternatives. For compounds 1 and 2, hydrogen bonding is only seen between the OH and NO groups; no hydrogen bonding is seen between OH groups. This limits the hydrogen bonding dimensionality to 0 or 1. It is found that compound 1 predominantly forms chains, while compound 2 forms chains in the enantiopure case and is capable of forming both dimers and chains in the racemic crystals. According to the Dreiding force field calculations, a chain hydrogen bonding network provides greater stability than a dimer hydrogen bonding pattern, which may partially explain the incorrect stability ranking of the crystal structure of compound 2. The aim of this study was to investigate whether crystal structure prediction approaches can be used to predict spontaneous resolution of chiral molecules. For the four compounds studied so far (the two compounds in this study and the two compounds in our previous study2), the simulations predict the correct stability order between enantiopure solids and racemic solids in three out of four cases. Advances in the accuracy of lattice energy calculations are required to improve this success rate. Acknowledgment. This project has been funded by the University of Bradford, School of Life Sciences. We thank Accelrys for providing a courtesy licence to the Materials Studio software package. Supporting Information Available: Crystallographic information files and four contour plots of energy versus T1 and T2 showing the results of conformational analyses of compound 1 and compound 2, along with complete lists of partial atomic charges for both compounds. This information is available free of charge via the Internet at http:// pubs.acs.org.

References (1) Collins, A. N.; Sheldrake, G. N.; Crosby, J. Chirality in Industry; John Wiley & Sons: Chichester, 1992; pp 5-54.

Crystal Growth & Design, Vol. 8, No. 8, 2008 2905 (2) Gourlay, M. D.; Kendrick, J.; Leusen, F. J. J. Cryst. Growth Des. 2007, 7, 56–63. (3) Ikuma, N.; Tamura, R.; Shimono, S.; Kawame, N.; Tamada, O.; Sakai, N.; Yamauchi, J.; Yamamoto, Y. MendeleeV Commun. 2003, 109– 111. (4) Ikuma, N.; Tamura, R.; Shimono, S.; Kawame, N.; Sakai, N.; Yamamoto, Y.; Yamauchi, J. Mol. Cryst. Liq. Cryst. 2005, 440, 23– 35. (5) Deumal, M.; Lafuente, P.; Mota, F.; Novoa, J. J. Synth. Meth. 2001, 122, 477–483. (6) Ikryannikova, L. N.; Ustynyuk, L. Y.; Tikhonov, A. N. J. Phys. Chem. A 2004, 108, 4759–4768. (7) Filippini, G.; Gavezzotti, A.; Novoa, J. J. Acta Crystallogr. B 1999, 55, 543–553. (8) Leusen, F. J. J.; Noordik, J. H.; Karfunkel, H. R. Tetrahedron 1993, 49, 5377–5396. (9) Leusen, F. J. J. Cryst. Growth Des. 2003, 3, 189–192. (10) Accelrys Cerius2, v 4.6; Accelrys Inc: San Diego, CA, 2001. (11) Accelrys Materials Studio, v 4.0; Accelrys Inc: San Diego, CA, 2006. (12) Guest, M. F.; Bush, I. J.; Van Dam, H. J. J.; Sherwood, P.; Thomas, J. M. H.; Van Lenthe, J. H.; Havenith, R. W. A.; Kendrick, J. Mol. Phys. 2005, 103, 719–747. (13) Allen, F. H. Acta Crystallogr. B 2002, 58, 380–388. (14) Karasawa, N.; Goddard, W. A. J. Phys. Chem. 1989, 93, 7320–7327. (15) Williams, D. E. Acta Crystallogr. A 1971, 27, 452–454. (16) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897–8909. (17) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (18) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785– 789. (19) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654– 3665. (20) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650–654. (21) Baur, W. H.; Kassner, D. Acta Crystallogr. B 1992, 48, 356–369. (22) Naik, N.; Braslau, R. Tetrahedron 1998, 54, 667–696. (23) Forrester, A. R.; Thomson, R. H. Nature 1964, 203, 74–75. (24) Laio, A.; Gervasio, F. L.; VandeVondele, J.; Sulpizi, M.; Rothlisberger, U. J. Phys. Chem. B 2004, 108, 7963–7968. (25) Neumann, M. A.; Leusen, F. J. J.; Kendrick, J. Angew. Chem., Int. Ed. 2008, 47, 2427–2430.

CG701256E