Effect of Structurally Similar Additives on Crystal Habit of Organic

Article pubs.acs.org/crystal

Effect of Structurally Similar Additives on Crystal Habit of Organic Molecular Crystals at Low Supersaturation Zubin B. Kuvadia and Michael F. Doherty* Department of Chemical Engineering, University of California Santa Barbara, California 93106-5080, United States ABSTRACT: Foreign molecules such as additives or impurities may influence crystal morphology to a significant extent by disrupting the growth mechanism and inhibiting the growth rate of certain crystal faces. The additives can be byproducts of a preceding reaction step, or they could be specially designed to obtain specific crystal morphologies for particular applications. Because of the specific chemical structure of structurally similar additives, they can incorporate or lock into the original lattice only in certain configurations; hence, only some of the crystal faces are able to recognize these additive molecules. Once attached, they interfere with the first turn of a rotating growth spiral and hence slow down the growth rate of the face. In this article, we develop a generic probabilistic scheme for quantitatively estimating imposter recognition on each crystal face using a combination of mean field theory and configurational energy minimization. On the basis of this recognition, the mechanistic effect of an impurity on the first turn of the spiral and hence the modified growth rates and modified crystal habits are computed. These concepts are generalized for all molecular crystals, including non-centrosymmetric molecules. We demonstrate the applicability of the model by correctly predicting the experimental morphologies of α-glycine with L-alanine impurity and paracetamol with p-acetoxyacetanilide impurity grown from water.

1. INTRODUCTION AND MOTIVATION Crystallization from solution is an essential unit operation in a wide range of chemical processing industries, e.g., catalysts, pharmaceutical compounds, food processing, man-made fibers, specialty and bulk chemicals. The crystallization process and operating parameters govern many fundamental properties of the resultant downstream materials including chemical purity profile, polymorphic state, crystal size and shape distributions, etc. Crystal shape or crystal morphology impacts the end-use efficacy of solid products (e.g., bioavailability for pharmaceutical compounds,1 reactivity for catalysts2). As the morphology of crystallized particles broadly determines their flow properties, appropriate morphology is essential for blending and compaction of drug particles. Filtration and drying times are also a function of crystal morphology; hence the economics of the downstream processing depends on the crystal morphology. Foreign molecules can be present in the solution from which crystals are being grown. Sometimes, foreign molecules such as impurities are unavoidable due to chemical pathways and selectivity constraints of the preceding reaction unit operation, whereas, in other cases, tailor-made additives are used to achieve specific objectives. Doping is also very common in the materials industry to achieve specific material properties. It is well-known from experiments that foreign molecules have a significant influence on the surface structure as well as morphology of crystals. There are numerous experimental efforts to © 2013 American Chemical Society

explain the effect of additives on crystal morphology of the native compound. Some of the most notable recent work in this field includes the study of the effect of Mg2+ ions on calcite by De Yoreo and co-workers3 or the work on impurity design and recognition for barium sulfate by Davey and co-workers.4 It was always observed that additives selectively affect certain faces of the crystal based on their chemical structure. This property was also used for separation of chiral compounds such as 5,6 L-alanine and D-alanine, L-leucine and D-leucine or to design tailor-made additives to achieve a particular morphology, for example, changing the morphology of the undesired L-cystine kidney stones.7 Previous modeling efforts in addressing the effects of additives on crystal morphology fall mainly into two broad categories: nonmechanistic and mechanistic. Nonmechanistic models such as the modified attachment energy approach8,9 are useful in detecting which faces are capable of recognizing an impurity, but give no mechanistic insights on how the growth rate of the face is inhibited. Kink density reduction models of Chernov,10 and Cabrera and Coleman,11 the step-pinning model of Cabrera and Vermilyea12 and the “spiral pinning model” of Sizemore and Doherty13 are the main mechanistic models. The kink density reduction model reasons that some of Received: July 25, 2012 Revised: February 2, 2013 Published: February 8, 2013 1412

dx.doi.org/10.1021/cg3010618 | Cryst. Growth Des. 2013, 13, 1412−1428

Crystal Growth & Design

Article

the kink sites are killed because of additive locking in those sites; however, it does not explain the effect of a supersaturation dead zone (i.e., the stop-start behavior of growth rate) seen at low additive concentrations. The step-pinning model proposed that the immobile impurities partition the edge into multiple segments, arresting the growth of those segments whose length is less than or equal to some critical length, thus decreasing the overall velocity of that edge. Some experimental papers14 contradict the step-pinning model. In this work, we build on and extend the spiral pinning model introduced by Sizemore and Doherty (2009).13 Our main focus is on molecular organic systems and their crystal shape evolution in the presence of structurally similar additives. A structurally similar additive is a foreign molecule present in solution that has many of the same atoms in much the same chemical structure as the host molecule. Such an additive usually differs from the host molecule by only one or two functional groups. Structurally similar additives are also called “tailor made additives” or “imposters” in the literature. The terms “structurally similar additive”, “additive”, “imposter” and “impurity” are used interchangeably in this work and are taken to mean the same thing in the context of crystal growth. The main goal of our method is to develop a rapid estimate to screen additive candidates (in conjunction with experiments) to engineer crystal morphology. The method is designed to give results in minutes rather than days. To achieve this we employ a mechanistic model which requires only very modest all-atom molecular mechanics and energy minimization calculations. The theoretical concepts used to develop these models rely on a mechanistic understanding of relevant surface phenomena governing the crystal growth process as well as knowledge of bulk solid-state chemistry. In particular, we develop a quantitative impurity recognition method on each crystal face using a combination of mean field theory and configurational energy minimization. On the basis of this recognition, we compute modified crystal growth rates of all facets and hence modified steady-state crystal shapes. We describe the challenges for non-centrosymmetric growth units and the methodology to overcome these challenges and generalize the methodology to all molecular crystals. We demonstrate the method for two systems: glycine/L-alanine and paracetamol/p-acetoxyacetanilide (PAA) grown from solution.

vi = 0 vi = v∞ , i

h τ

N

∑ i=1

3. IMPURITY RECOGNITION MODEL As stated in the Introduction, it is experimentally observed that only certain crystal faces are affected by the impurity/additive (this changes the relative growth rates of all faces and thus modifies the crystal shape); hence we need a quantitative method to estimate the composition of the impurity that is incorporated on each crystal face. The term incorporation in this context means that the impurity/additive is locked into the solid in a host kink site along a step and is included as an impurity in the solid crystalline phase. This composition should depend on the composition of impurity present in the solution from which the crystals are growing as well as the energetics of the crystal system and the solution. At low supersaturations, where spiral growth is the dominant mechanism, Sizemore18 makes a mean-field lattice approximation to estimate the free energies of the terraces of different orientations hkl and the solution. Assuming the additive to be distributed in equilibrium between the solution and the crystal terrace, an equation was proposed to estimate the incorporation composition (mole fraction) of additive on that terrace, zp, by equating the chemical potentials of the additive on the terrace and in the solution 2⎫ ⎧ ⎪ 2/ vp(l) − / vv(l) − (2/hp − /hh)(1 − z p) ⎪ ⎬ z p ≈ xp exp⎨ ⎪ ⎪ kBT ⎭ ⎩

(5)

where xp is the mole fraction of the imposter in solution, / vp(l) represents the total interaction energy between a molecule of imposter and a sea of solvent, / vv(l) represents the total interaction energy between a solvent molecule and a sea of itself, / hp represents the total interaction energy between an imposter and a terrace of solute host molecules of orientation (hkl), and / hh represents the total interaction energy between a solute host molecule and a terrace of itself of orientation (hkl). The symbol l signifies “liquid” whereas v, p, and h signify solvent, impurity and host (solute), respectively. The algorithm to compute all the interaction energies above is described in section 3.1. There can be multiple types of growth units within a single unit cell (which may either be related to each other by symmetry operators or not), each one of which can be potentially replaced by the structurally similar additive, which brings up the question of selecting the most favorable positions for the additive molecule on a given slice.

(1)

lc , i + 1 sin(αi , i + 1) vi

(3)

On the basis of eqs 1 and 2, a decrease in growth rate of a crystal face can be caused either by a reduction in step velocity (as proposed by the Cabrera and Vermilyea step pinning mechanism12 and the kink density reduction mechanism10,11) or by an increase in the critical length of edges on that face (discussed in the spiral pinning mechanism proposed by Sizemore and Doherty13).

where τ is the characteristic rotation time of the growth spiral, which corresponds to the time between consecutive passes of steps across a point on the face, and h is the step height.15 Snyder and Doherty15 give an expression for the rotation time of a convex N sided polygonal spiral as τ=

li > lc , i

whereby a spiral edge starts growing with a constant step velocity only after reaching a critical length. Since kink sites are the incorporation sites for solute molecules, the step velocity is the product of the distance of propagation (ap,i), density of the kink sites (ρi), and the rate of incorporation into the kink sites (ui). Hence, vi ∝ ap , iuiρi (4)

2. CRYSTAL MORPHOLOGY MODEL AT LOW SUPERSATURATION The growth rate of a crystal face that is growing by the spiral mechanism can be expressed as G=

li ≤ lc , i

(2)

where αi,i+1 is the angle between edges i and i+1; lc,i+1 is the critical length of edge i+1; and vi is the velocity of step i. We use the step velocity profile first proposed by Kaischew and Budevski16 and later by Voronkov17 1413



Article

Figure 1. Diagram illustrating native and modified slice and attachment energies. Reproduced from Clydesdale et al.8 with permission. Copyright 1994 Elsevier.

where n represents the total number of different growth units within a single unit cell comprising a single surface layer that can be replaced by an imposter molecule. (All positions within a unit cell may or may not be present on a given crystal face hkl based on extinction symmetry conditions; faces in which the interplanar spacing is half the unit cell distance in that direction will consist of only half the total number of growth units of the unit cell. For example, on the {020} face of α-glycine, only two docking positions exist, whereas there are four positions in the unit cell and on other faces of this crystal. For each face we select the positions possible for a given face.) This approach ensures that the probability of docking in each of the most favorable positions (aligned positions) on a face is correctly accounted for. Note that each of the mole fractions zp,i in eq 6 are calculated using eq 5 with the appropriate value of / hp for the imposter placed in the ith position in the unit cell. Note that / hp conceptually differs from delta binding energy (Δb), while the latter signifies the compatibility of an additive to the lattice position and its likelihood of incorporation, / hp is a measure of the distortion of the lattice caused by the additive once it has incorporated into the lattice. The quantitative difference between the two is pointed out in section 3.1 3.1. Algorithm To Find the Interaction Energies of Solid Terrace and Solution Phase. In order to compute the delta binding energies of different additive positions, we need the modified slice and the modified attachment energies when each growth unit in the unit cell is replaced in turn by the additive growth unit as well as the original host−host lattice interaction energies in the absence of any additive. The solid− solid intermolecular interactions are computed using the program HABIT19 and applying the AMBER force field.20,21 The partial atomic charges on each of the atoms of all the molecules comprising the unit cell were calculated using GAUSSIAN 0322 with the restrained electrostatic potential (RESP) model.23 The quantity / hh is calculated by the sum of the solid−solid bonds

In principle, the additive can dock in any position (with a corresponding probability) in a given slice. However, we make the simplifying assumption that the most favorable additive positions are those which are completely aligned to the original host molecules. Our method is consistent with previous methods in the field of crystal growth. The magnitude of the deformity of the original lattice structure when one host molecule is replaced by the additive is given by the delta binding energy (Δb) concept.8,9 In other words, the delta binding energy is used to assess how easily an additive will incorporate into a given crystal surface and is a measure of the molecular compatibility of the host and additive system. The delta binding energy is the difference between the incorporation energy of the additive and the pure host. The delta binding energy compares the energies of a host molecule after being bound up in an (hkl) layer (Esl) and subsequent attachment of the layer to a substrate (Eatt) to those of an imposter molecule bound up in the same (hkl) layer (Esl′ ) and subsequent attachment of that af fected layer to a substrate ′ ). Defining Δb ≡ (Esl′ + Eatt ′ ) − (Esl + Eatt), the lower the (Eatt value of Δb, the more favorable the additive binding. (We use the attachment energy definition of Berkovitch-Yellin et al.9 who define it as the energy per molecule released when a layer is attached to the crystal lattice; see Figure 1 for illustration of different terms). The net mole fraction of the additive on the surface is taken as the weighted average of the mole fractions found at each of the individual positions, zp,i, using the Boltzmann factor of the binding energy (i.e., binding probabilities) for averaging, given by the following equation exp zp =

−Δb1 kT

−Δb2 kT

−Δbn kT

( )z + exp( )z + ......exp( )z exp( ) + exp( ) + ......exp( ) p,1

−Δb1 kT

−Δb2 kT

p,2

p, n

−Δbn kT

(6) 1414



Article

(3) This structure is still not the optimal docking (aligned) position for the additive, and for most molecules we must use restrained configurational energy minimization to get the optimal structure of additive (i.e., the optimal position of all its atoms) that can be used for solid-state energy calculations. We replace one of the host molecules on any chosen hkl surface (F-face) with the additive molecule, as determined by steps (1) and (2). The imposter is not restrained and has all degrees of freedom (to relax any deformed bonds or angles), whereas all the other host molecules on the surface and the bulk are position restrained (since they are expected to sit in their already optimized crystal lattice geometrical positions). We use AMBER10 package20 and atomistic simulations with a large number of cycles (100 000) to carry out energy minimization (7000 cycles of steepest gradient, the rest using conjugate gradient). We used a cutoff based method (no periodicity is applied and Particle Mesh Ewald was off; a cutoff of 1.5 nm was used. The molecular system considered to represent the particular face was 1085 atoms, 53 host molecules and 1 additive molecule). The crystal face was explicitly generated by feeding in geometrical coordinates from the crystal visualization software Mercury. Energy minimization was carried out for just one position on the chosen slice. The end configuration for the additive molecule (with the application of appropriate symmetry operators to replicate all the n different favorable positions in the unit cell for direct substitution in different slices) can be used for computing the energetics (/ hp) to enable estimation of the net composition (mole fraction) of additive on each face hkl of the crystal. The solution interaction energies are estimated using bulk thermodynamic quantities such as enthalpy of vaporization for / vv(l) and the enthalpy of solvation for / vp(l). More details are available in Appendix B of Sizemore(2008).18 We can now directly use eqs 5 and 6 to find the net concentration of additive on any given face hkl.

satisfied upon forming a slice of orientation (hkl) and the solid− solid bonds satisfied upon attaching that slice to an existing surface of orientation (hkl). The procedure is similar for the additive for the calculation of / hp, whereby the central molecule is replaced by the additive molecule with optimized positional coordinates for each of the n positions it occupies in the host crystal lattice. The modified slice (Esl′ ) and the modified ′′ ) are computed for each of the replaced attachment energies (Eatt positions. Hence, we get /hh = Esl + Eatt ″ /hp = Esl′ + Eatt

(7)

(Note: As mentioned in the previous section, the quantity / hp is different from the incorporation energy used for the delta binding energy (Δb) concept (Esl′ + Eatt ′ ). The quantity (E′sl + E′att′ ) is the measure of the change in the lattice energy caused by the additive molecule af ter it has docked in a given position and additive is the central molecule for the energy ′ ) is the change caused because calculation, whereas (Esl′ + Eatt of a slice containing the additive docking into a pure lattice in which the host is the central molecule, see Figure 1 for illustration of E′att and E′att′ ). We need to compute the best alignment of the additive with the host molecule in order to substitute additive molecules for host molecules in the lattice. This is an important step as it forms the initial condition for computing all the modified solid−solid interactions. We can start off directly with any available Cartesian coordinates of the additive molecule; however we chose the unit cell crystal structure of the additive molecule as a convenient set of starting coordinates for the additive molecule. It is necessary that the volume of the additive growth unit is preserved while transforming its coordinates to the host lattice system; i.e., the intramolecular bonds and bond angles are not deformed. The overall scheme that we use is comprised of three steps: (1) The first step is to find the optimal rotation matrix using the Cartesian coordinates of the common atoms. The common atoms are the atoms of the fragment that the host molecule and the structurally similar additive have in common with each other. The structurally similar additive usually differs from the host molecule by the presence of a single functional group. Several fast optimal superposition (FOS) techniques are available24,25 and are used in the field of computational structural biology, from which we borrow the general methodology which is based on linear algebraic operations. We compute the singular value decomposition of the correlation matrix (which is generated using the identical (common) atoms of the host and the additive) adapting the method described by Lesk (1986).26 The rotoinversion, if detected while finding the optimal superposition, is eliminated by changing the sign of the singular vector corresponding to the minimal singular value. The singular value decomposition directly gives us the optimal rotation matrix while simultaneously preserving the volume of the additive molecule (verified by the determinant of the matrix being equal to +1). In this operation, the additive (the common fragment of the additive) is considered to be rigid; hence we are simply rotating the whole molecule to ensure the volume preservation. (2) We then minimize the RMSD of the common atoms between the two aligned molecules to find the optimal translation matrix to superimpose the host and the additive so that the additive can replace one of the host molecules in the host lattice.

4. MODELING THE EFFECT OF ADDITIVES ON FIRST TURN OF SPIRAL According to the spiral pinning model, additives have no effect on fully developed spiral edges, since they have already reached their critical lengths and move with a constant velocity. The effect is seen only on the first turn of an emerging spiral where the sides are yet to reach their critical lengths. The additive, if present on that particular face, increases the rotation time of the spiral by increasing the critical lengths of all edges, i.e., causing the edges to take longer to reach constant normal velocity. Thus, it reduces the growth rate of that particular face. The total time required for edge i to elongate such that it can begin its own growth is13 τi =

z p(1 + mci ) ⎤ lc, i sin(αi − 1, i) ⎡ ⎢1 + ⎥ , z p ∈ [0, z p*) 1 − z p(1 + mci ) ⎥⎦ vi − 1 ⎢⎣ (8)

where mci ≡ lc,i/aei, which is the critical length expressed in multiples of the lattice parameter along edge i. Typical values for mci are on the order of 100. To calculate the critical lengths, we use the Gibbs−Thomson Law applied to a single spiral edge.15,27,28 This law is a thermodynamic condition which states that adding a row of molecules of length li greater than the critical length lc decreases the free energy of the system, whereas adding a row less than the critical length increases the free energy of the system; thereby 1415



Article

preventing its growth. Thus, the critical length is the length at which the change in free energy of adding a new row of molecules to the system is zero, and can be calculated from the kink energy of the given edge using the equation28

lc, i =

Table 1. Relative Growth Rates of Crystal Faces of α-Glycine Grown from Aqueous Solution: No Additive Presenta ρ

lc (Å)

dhkl (Å)

R

(020)

[001] [100] [001] [100] [100] [111̅] [001] [111̅]

0.035 0.162 0.035 0.162 0.162 0.181 0.035 0.181

895 489 895 489 489 668 895 668

5.98

1.26

5.98

1.26

4.67

3.38

4.41

1

(020̅ ) {011} {110}

(9)

2ae, i⟨ϕkink, i⟩ RT ln S

(10)

where Δμ is the difference between the chemical potential of the solute molecules in solution and in the crystal, S is the supersaturation ratio, ae,i is the distance between molecules (centers of mass) along edge i and ⟨ϕkink,i⟩ is the arithmetic average kink energy on edge i. The kink energy can be computed using the classical BCF method as the energy required to form kinks from a straight step by the process of rearrangement.15,29 We see that the critical length has inverse dependence on supersaturation. This means that at high supersaturation, lc,i has decreased such that a given zp is less effective at slowing growth. Hence, a higher supersaturation diminishes the effect of additives, whereas at very low supersaturation, the effect of additives is more pronounced (resulting in a dead zone of zero growth of that face).

Figure 2. (a−c) Steady-state morphology of α-glycine (three views).

bond chains

Δμ

which then becomes lc, i =

face

2ae, i⟨ϕkink, i⟩

5. APPLICATION TO SYSTEMS OF REAL-COMPLEXITY 5.1. α-Glycine plus Additive L-Alanine Grown from Aqueous Solution. The system of α-glycine with L-alanine was specifically selected to demonstrate the anisotropic or the polar effect an additive can have on the crystal morphology. The α polymorph of glycine (zwitterionic) crystallizes from water in the space group P21/n with unit cell parameters given by a = 5.105 Å, b = 11.969 Å, c = 5.465 Å, and β = 111.70°; there are four molecules in the unit cell, with Z = 4 symmetrically equivalent positions and J = 1 molecule per equivalent position. On the basis of the literature which suggests that glycine exists as dimers in solution30 and on AFM

a

The {110} faces are taken to be the reference for calculating relative growth rates.

Figure 3. Four different positions for the imposter L-alanine to replace glycine molecules and form dimers within the unit cell (scales are in Å). 1416



Article

Figure 4. The net mole fraction of the additive (L-alanine) on different slices of α−glycine with varying mole fractions of the additive in solution. Only the (020̅ ) face has any detectable incorporation on the surface. The impurity does not affect (is not recognized by) any other face.

Table 2. Additive Mole Fractions in Solution, xp, and on the (020̅ ) Face, zp, as a Function of Supersaturation

studies which indicate that step heights on the crystal terraces are equivalent to dimers,31 we assume that the glycine molecules incorporate and disincorporate from the crystal lattice as dimers (Z = 2, J = 1). The dimer growth unit is centrosymmetric, and hence the kink rate, ui, in eq 4 is the same for every edge and can be omitted. We predict the steady-state morphology of glycine from water to be the familiar coffinshape as shown in Figure 2, which is in good agreement with the famous experimental shape by Weissbuch et al.5,6 The relative growth rate results for α-glycine are presented in Table 1 (R represents the relative growth rate of face hkl normalized by the reference {110} faces; ρ and dhkl signify the kink density for the edge and the interplanar spacing of the face, respectively). Our predictions are made at a low supersaturation (4%). The solid state force field chosen was the Generalized Amber Force Field21 using RESP partial charges obtained from Gaussian03,22 which gave a calculated lattice energy of −68.0 kcal/mol. This compares favorably to the experimental lattice energy of this zwitterionic compound of −67.25 kcal/mol32 (computed using the experimental values of the sublimation enthalpy (34.7 kcal/mol33), proton transfer energy (31.34 kcal/ mol32) and adjusted by a factor of αRT where α is typically taken as 2.) All the solid−solid broken bonds used in the calculation of kink energies need to be interfaced to account for the effect of the solvent. We scale the bonds according to the approach described by Winn and Doherty34,35 and extended by Snyder and Doherty.15 The method is a classical one based on estimating solid and solvent interactions independently (making use of solubility parameters for the solvent phase36) and calculating the interfacial energy from the adhesive and cohesive components. The chiral nature of the additive calls for subtlety while using the singular value decomposition method to prepare the initial imposter structure for various symmetry-related positions in the lattice as certain symmetry operations (e.g., reflections, inversions) cannot be carried out for a chiral compound.

S

conc of glycine (g/L)

conc of L-alanine (g/L)

xp

zp

1.05 1.1 1.15

262.5 275 287.5

2.625 2.75 2.875

0.000499 0.000521 0.000544

0.0070 0.0071 0.0072

Table 3. Modified Growth Rates of Faces of Glycine in the Presence of 1 wt % L-Alaninea face

S

zp

G(zp)

G(0)

G(zp)/G(0)

R(zp)

(020)

1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15

0 0 0 0.0070 0.0071 0.0072 0 0 0 0 0 0

0.000821 0.00164 0.002461 0 0.00091 0.00173 0.0022 0.00439 0.00659 0.00065 0.001306 0.00195

0.000821 0.00164 0.002461 0.000821 0.00164 0.002461 0.0022 0.00439 0.00659 0.00065 0.001306 0.00195

1 1 1 0 0.555 0.703 1 1 1 1 1 1

1.26 1.26 1.26 0.00 0.70 0.89 3.38 3.38 3.38 1.00 1.00 1.00

(020̅ )

{011}

{110}

a

The {110} faces are taken to be the reference for calculating relative growth rates.

In symmetry positions that use the reflection or the inverse operation to generate molecules, we create D-alanine using the optimal superposition matrix before applying the reflection symmetry operator to end up with L-alanine molecules at all symmetry positions in the lattice after their respective symmetry operations have been applied (see Sizemore (2008)18 pp 109− 110 for additional details). A single molecule of L-alanine has an extra methyl functional group compared to glycine and is thus 1417



Article

As described in section 3.1, we first align the L-alanine additive molecule to the host glycine molecule and then replace each host growth unit in the unit cell sequentially with the corresponding additive growth unit. We calculate the modified slice and attachment energies of all possible positions on a given face and then compute the binding energies. The binding energy values are shown in Table 7, and the solid state terrace interaction energies for the host and the additive molecule in different positions on different faces are shown in Table 8. As mentioned earlier, the {020} face contains only two docking positions, whereas there are four positions in the unit cell and on other faces. Knowing the composition of the additive in the solution, we compute the net mole fraction of the additive on each face using the mean field theory (eq 5) and the binding probability (Boltzmann factor of binding energy) for each of the n possible positions on the face hkl (eq 6). The solution interaction energies for this system are taken from the experimental values reported for water (/ vv(l) for water is −9.7 kcal/mol whereas / vp(l) for the glycine/L-alanine additive dimer growth unit is −42.5 kcal/mol37). Figure 4 shows the additive mole fraction on different faces of α-glycine as a function of the composition of the additive in solution. Only the (020̅ ) face is influenced by the additive, and all the other faces (including the opposite face (020)) remain totally unaffected by the additive. In other words, the other faces do not recognize the L-alanine additive, whereas the additive mole fraction on the (020̅ ) face steadily increases with increasing solution mole fraction of the additive.

Figure 5. The predicted crystal morphology of α-glycine grown from aqueous solution in the presence of L-alanine at different supersaturation, at a constant temperature of 25 °C. Additive composition is 1% by mass of the solute in each case, which corresponds to different imposter mole fractions in solution as shown in Table 2.

capable of forming hydrogen bonds with a glycine host molecule in solution (strong hydrogen bonds form between carboxylic and amino groups), so we assume that the glycine/L-alanine dimer is the crystallographic growth unit that attaches and detaches from kink sites. There are four different potential docking sites for L-alanine within the α-glycine unit cell that comprises two dimer units of α-glycine (see Figure 3).

Figure 6. Experimentally grown α-glycine crystals in aqueous solution in the presence of L-alanine (1 wt % of solute) at 28 °C and ≈5% supersaturation (solubility of glycine at 28 °C = 266 g/L; glycine concentration in solution = 280.6 g/L; L-alanine concentration in solution =2.85 g/L; xp = 0.000539). (a) and (b) are optical images, whereas (c) is an image obtained using SEM. The arrows point to crystals with the classical half coffin shapes seen among a dispersion of crystal shapes. The half coffin shape is also clearly visible in (c). 1418



Article

With the calculated values of zp for L-alanine, the corresponding imposter-modified growth rates for α-glycine and hence its modified crystal morphology can be computed. We show the calculations for three different supersaturation values at a constant temperature of 25 °C. In each case, the additive is 1 wt % of the different solute concentrations (see Tables 2 and 3, the solubility of α-glycine in water at 25 °C is 250 g/L). As seen in Table 3, the (020̅ ) face exhibits a dead zone at low supersaturation, reflected by zero growth rate at S = 1.05. This is because this mole fraction of additive on the face prevents at least one edge from ever reaching its critical length. Hence, this edge cannot start growing, thus pinning the spiral growth in its first turn. The dead zone is released at higher supersaturation because the critical length decreases with increase in supersaturation, so now the edges are able to grow out and complete the spiral turns. Another interesting point to note from Table 3 is that the relative growth rates of faces in the absence of any additive are constant with increasing supersaturation, as anticipated in the spiral growth regime.38 The additive composition on the terrace is almost constant at all three supersaturations. This behavior is discussed in section 6. The predicted crystal morphologies in the presence of additive at the three different supersaturations are shown in Figure 5. Note that we repeated all these calculations using α-glycine monomer as the growth unit and found that the additive is again recognized exclusively on the (020̅ ) face. Thus, our model agrees with the experimental observations regardless of which type of growth unit is selected. It is well documented from experiments6 that the chiral additive, L-alanine, affects only the (02̅0) face of the native α-glycine crystal and not its opposite (020) face which belongs to the same family. We performed our own experiments of growing α-glycine crystals in the presence of L-alanine from aqueous solution by the cooling crystallization method. The crystals were grown at low supersaturation in the range 5−8%, at 28 °C in a small, quiescent crystallizer equipped with a video microscope, moveable stage, and Peltier cell for precise temperature control (experimental setup is described elsewhere39,40). The composition of the additive in solution was 1% by mass of the solute, the same as with the predictions. Figure 6 shows the experimental morphologies: (a) and (b) are the optical microscope images, whereas (c) is an image captured using a scanning electron microscope (SEM). The arrows point to crystals with the classical half-coffin shape that results because of the additive. There is a dispersion in crystal shape which is common during experiments because of local variations in supersaturation due to diffusion effects and also due to fluctuations in the composition of additive uptake of each crystal. It is difficult to maintain constant and uniform supersaturation during crystallization experiments as discussed in section 6. 5.2. Paracetamol plus PAA Grown from Aqueous Solution. Paracetamol (acetaminophen) is an important and widely used analgesic drug. Paracetamol has been widely studied in the literature because of its significance in the pharmaceutical industry and because of its well-faceted crystal habit in different solvents. There have been various attempts to study the shape evolution and develop models to predict its crystal shape. Monoclinic paracetamol crystallizes in the space group P21/a. There are four molecules per unit cell related to each other by symmetry operations; the cell parameters are a = 12.92 Å, b = 9.40 Å, c = 7.10 Å and β = 115.9°. Paracetamol, being a noncentrosymmetric molecule, has its own set of challenges in predicting crystal morphology, either in the absence or presence of impurity. Hence, this system is chosen

Figure 7. The predicted steady-state morphology of paracetamol grown from aqueous solution in the absence of PAA using our model for noncentrosymmetric growth units at low supersaturation (4%). Adapted from Kuvadia and Doherty38 with permission. Copyright 2011 America Chemical Society.

Figure 8. The starting structures of paracetamol (left) and PAA (right) taken from their respective crystal structures viewed in the same direction (perpendicular to the b axis). The paracetamol and PAA molecules differ by an acetoxy functional group (highlighted). C = carbon, H = hydrogen, O = oxygen, N = nitrogen.

Figure 9. Alignment of PAA (right) with paracetamol (left) after optimal rotation. Note the change in atomic numbering of PAA from its starting structure to match that of paracetamol as we now create an imposter structure that is compatible with the host lattice.

as the second case study to test the applicability of our model to complicated growth units. The challenges for this class of molecules for the morphology prediction in the absence of additives was discussed by Kuvadia and Doherty (2011).38 In that work, the crystal morphology of paracetamol was predicted using the spiral growth model for non-centrosymmmetric 1419



Article

Figure 10. One of the docking positions of PAA in the {110} face after configurational energy minimization of the additive molecule. The other positions on this face and other faces are related by symmetry to this position. Note: this figure shows only part of the simulation box. The simulation procedure is described in 3.1, step 3.

Table 4. Reduced Growth Rates of Different Faces of Paracetamol in the Presence of PAA at Three Different Supersaturation Values at a Constant Temperature of 28 °Ca face 110

001

011

201̅

020

200

S 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15

G(0) 0.000006753 0.000027070 0.000061110 0.000011830 0.000047330 0.000106590 0.000014218 0.000056890 0.000128110 0.000018270 0.000072980 0.000164000 0.000378800 0.001508000 0.003300000 0.000167900 0.000671600 0.001511000

R(0) 1.00 1.00 1.00 1.75 1.75 1.74 2.11 2.10 2.10 2.71 2.70 2.68 56.09 55.71 54.00 24.86 24.81 24.73

xp 3.57477 3.74463 3.91446 3.57477 3.74463 3.91446 3.57477 3.74463 3.91446 3.57477 3.74463 3.91446 3.57477 3.74463 3.91446 3.57477 3.74463 3.91446

× × × × × × × × × × × × × × × × × ×

−5

10 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−5

zp

G(zp)

G(zp)/G(0)

R(zp)

0.0035 0.0036 0.0037 0 0 0 0.000054 0.000057 0.00006 0 0 0 0 0 0 0.000034 0.000036 0.000038

0.000001944 0.000017300 0.000045110 0.000011830 0.000047330 0.000106590 0.000014040 0.000056480 0.000127700 0.000018270 0.000072980 0.000164000 0.000378800 0.001508000 0.003300000 0.000166600 0.000670000 0.001509000

0.288 0.639 0.738 1.000 1.000 1.000 0.987 0.993 0.997 1.000 1.000 1.000 1.000 1.000 1.000 0.992 0.998 0.999

1.00 1.00 1.00 6.09 2.74 2.36 7.22 3.26 2.83 9.40 4.22 3.64 194.86 87.17 73.15 85.70 38.73 33.45

a

The additive is 2 wt % of the solute in each case. The relative growth rates are defined with respect to the {110} faces. R(0) and R(zp) are relative growth rates in the absence of and in the presence of additive, respectively.

organic molecules which is shown in Figure 7 and is in good agreement with the experimental shapes reported by Ristic and co-workers.41 PAA has been used as a synthetic impurity for this drug to study the effect of a structurally similar additive on various crystal properties.42 Figure 8 shows the starting structures of paracetamol and PAA. We use the methodology described in the above sections to estimate the impurity recognition on different faces and its effect

on growth rates of faces and on crystal morphology. As there are four molecules within the unit cell related by symmetry operators, there are four possible sites within the unit cell where the additive PAA molecule can dock. The first step is to find the best alignment of an additive PAA molecule with one of the host paracetamol molecules. This is done in three steps as described in section 3.1. (1) Finding the optimal rotation using singular value decomposition (see Figure 9). 1420



Article

Figure 11. The predicted steady-state morphology of paracetamol grown from aqueous solution using our model for non-centrosymmetric growth units in the presence of 2% PAA at low supersaturation (5%).

Figure 12. The predicted habit evolution of a paracetamol seed grown from aqueous solution using our model for non-centrosymmetric growth units in the presence of 2 wt % PAA at low supersaturation (5%). This figure shows the shape evolution with the dimensionless time for four faces (a) {001}, (b) {201̅}, (c) {011}, (d) {110}. x represents the dimensionless distance of that face from the center of the crystal whereas the independent variable is the number of time steps reached in the integration of the differential shape evolution equation, Δξ = 0.25, ξmax = 25 (see Zhang et al.46 for details of the shape evolution equations). The graph in (a) is equivalent to the aspect ratio of the crystal.

different slices. We calculate the modified slice and attachment energies of all the possible positions and also compute the binding energies using the program HABIT19 and applying the AMBER Force field.20,21 The partial atomic charges on each of the atoms of the four molecules comprising the unit cell needed for the force field were calculated using GAUSSIAN 0322

(2) The optimal translation by minimizing RMSD to superimpose the molecules. (3) Configurational energy minimization to find the best alignment of the additive molecule in the host lattice (see Figure 10). We use this starting configuration of the additive molecule and use appropriate symmetry operators to replicate it in 1421



Article

Figure 13. Experimental shapes of paracetamol crystals grown in aqueous solution at T = 28 °C with 2 wt % PAA and a superaturation of 5%. (The saturation concentration of paracetamol at this temperature is 16.11 g/L.) The mole fraction of PAA at this conditions is 3.6 × 10−5. (a) and (b) are optical images, whereas (c) is an image obtained using SEM.

The calculations for a specific constant temperature of 28 °C at three different supersaturations are shown in Table 4. The saturation concentration of paracetamol at this temperature is 16.11 g/L.44 The additive is taken to be 2 wt% of paracetamol in each case. The solution mole fractions, xp, at each of three supersaturations are also presented in Table 4. The {110} faces are the most affected by the imposter followed by {011} and {200} faces. These are the faces in which PAA can dock in some of the lattice positions and act as an imposter with its acetoxy functional group projecting out of the slice (see Figure 10). The predicted morphology of paracetamol in the presence of 2 wt% PAA (of paracetamol) at 5% supersaturation is shown in Figure 11. The influence of the additive slows down the growth rate of the {110} faces significantly and hence the crystal habit is in the form of long bars. It can also be observed from Figures 11 and 12 and Table 4 that the predicted aspect ratio is about 6. (The aspect ratio of the native crystal in absence of additive is 2). As shown in Figure 12, the {201}̅ and the {011} faces, despite being F faces with strong periodic bond chains, disappear during the growth of the seed crystal and do not appear in the final steady-state predicted morphology.

with the restrained electrostatic potential (RESP) model,23 results are shown in Table 9. (Note: During the calculations, it is important to calculate the average slice and the attachment energies for the entire unit cell for both the host and additive calculations, rather than using the energy reported by HABIT for each individual position as described in detail in Appendix B of Kuvadia.43 This is significant for paracetamol because it is a non-centrosymmetric molecule and different host positions within the unit cell for a given slice will have different interaction energies.) Only those faces which have at least one lattice position with nonzero binding probability are shown in Table 9. Knowing the composition of the additive in solution, we use the mean field theory (eq 5) and the binding probabilities of all positions on these faces (eq 6) to compute the net mole fraction of the additive on each face. The solution interaction energies for this system are taken from the experimental values reported for water (/ vv(l) for water is −9.7 kcal/mol whereas / vp(l) for the PAA molecule is taken as −16.5 kcal/mol (using the ACD/PhysChemSuite software (Advanced Chemistry Development Inc.)). 1422



Article

Table 5. Reduced Growth Rates of Different Faces of Paracetamol in the Presence of PAA at Three Different Supersaturation Values at a Constant Temperature of 28 °Ca face

S

G(0)

R(0)

xp

zp

G(zp)

G(zp)/G(0)

110

1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15 1.05 1.1 1.15

0.000006753 0.000027070 0.000061110 0.000011830 0.000047330 0.000106590 0.000014218 0.000056890 0.000128110 0.000018270 0.000072980 0.000164000 0.000378800 0.001508000 0.003300000 0.000167900 0.000671600 0.001511000

1.00 1.00 1.00 1.75 1.75 1.74 2.11 2.10 2.10 2.71 2.70 2.68 56.09 55.71 54.00 24.86 24.81 24.73

0.00018 0.00019 0.00019 0.00018 0.00019 0.00019 0.00018 0.00019 0.00019 0.00018 0.00019 0.00019 0.00018 0.00019 0.00019 0.00018 0.00019 0.00019

0.0083 0.0085 0.0085 0 0 0 0.00026 0.00027 0.00027 0 0 0 0 0 0 0.00017 0.00018 0.00018

0.000000000 0.000000000 0.000027500 0.000011830 0.000047330 0.000106590 0.000013200 0.000054800 0.000125000 0.000018270 0.000072980 0.000164000 0.000378800 0.001508000 0.003300000 0.000161200 0.000657000 0.001488000

0.000 0.000 0.450 1.000 1.000 1.000 0.928 0.963 0.976 1.000 1.000 1.000 1.000 1.000 1.000 0.960 0.978 0.985

001

011

201̅

020

200

a

The additive is 10 wt% of the solute in each case.

Figure 14. The reduced growth rates of the (02̅0) face of α-glycine grown from aqueous solution in the presence of different L-alanine solution compositions at different supersaturations of glycine.

These faces do appear, however, in some experimental images, and we also see a dispersion in experimental crystal shapes (see Figure 13). We grow crystals of paracetamol in the presence of 2 wt% PAA (relative to the paracetamol) at low supersaturation (5−8%) and at 28 °C using the method described in section 5.1. Experimentally, we observe a range of aspect ratios (4−10) as seen in Figure 13. Hendriksen et al.45 report a mean aspect ratio of 6.7 (with a standard deviation of 4) at a similar PAA concentration (2.5 wt %) for paracetamol crystals grown by slow

evaporation (corresponds to low supersaturation). They also reported an aspect ratio of 6.84 (SD 2.58) for crystals grown by cooling. Hence, our predictions are in agreement with the experiments. We also predict an increase in aspect ratio of the crystal dimensions (length/width ratio, i.e., (001):(110)) with an increase in additive concentration and a decrease in aspect ratio with increase in supersaturation (as seen from Tables 4 and 5). This is in agreement with the experimental trends observed by Chow et al.42 Also, at higher additive composition, there is a 1423



Article

more pronounced dead zone reflected by zero growth rate at low supersaturation; see Table 5. The predicted aspect ratio should be very high in that case (long rods/needles), again in agreement with the trend seen in experiments.42

Table 6. The Mole Fraction of Additive, L-Alanine, (zp) on the (020̅ ) Face Vs Solution Mole Fractions, xp, at 5% Supersaturation L-alanine

wt% of glycine 1 2 4 6 10 15 20

xp

zp

0.0005 0.001 0.002 0.003 0.005 0.0075 0.0098

0.007 0.0088 0.0107 0.0118 0.0133 0.0146 0.0154

6. DISCUSSION This work advances the spiral pinning mechanism in which the additives slow down the rotation speed (increase the rotation time) of the first turn of a growing spiral. Our model gives a very reasonable quantitative estimate of which crystal facets are affected by a structurally similar additive and to what

Table 7. Binding Energies for Slices of α-Glycinea slice

pos

Esl

Eatt

E′sl

ΔEsl

ΔE′att

Δb

020̅

1 2 1 2 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

−78.5 −78.5 −78.5 −78.5 −36.264 −36.264 −35.823 -35.823 −35.823 −35.823 −36.265 −36.265 −35.823 −35.823 −36.265 −36.265 −36.264 −36.264 −35.823 −35.823 −50.611 −50.611 −50.728 −50.728 −50.728 −50.728 −50.611 −50.611 −50.728 −50.728 −50.611 −50.611 −50.611 −50.611 −50.728 −50.728 −17.784 −17.784 −17.784 −17.784 −17.784 −17.784 −17.784 −17.784

−6.865 −6.865 −6.865 −6.865 −24.725 −24.725 −24.725 -24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −33.855 −33.855 −33.855 −33.855 −33.855 −33.855 −33.855 −33.855

−52.85 −75.313 −52.85 −75.313 −32.834 518.63 −31.89 -35.579 −31.89 −35.579 −32.83 518.63 −31.89 −35.579 −32.83 518.63 −32.834 518.63 −31.89 −35.579 −39.43 393.48 −38.93 57.044 −38.93 57.045 −39.43 393.48 −38.93 57.045 −39.43 393.48 −39.43 393.48 −38.93 57.044 −18.38 90.565 −18.38 90.569 −18.38 90.565 −18.38 90.569

25.65 3.187 25.65 3.187 3.43 554.894 3.933 0.244 3.933 0.244 3.435 554.895 3.933 0.244 3.435 554.895 3.43 554.894 3.933 0.244 11.181 444.091 11.798 107.772 11.798 107.773 11.181 444.091 11.798 107.773 11.181 444.091 11.181 444.091 11.798 107.772 −0.596 108.349 −0.596 108.353 −0.596 108.349 −0.596 108.353

−0.34 −0.09 −0.17 447.06 −0.76 −1.54 12.78 -0.81 12.78 −0.81 −0.76 −1.54 −0.76 −1.54 12.78 −0.81 12.78 −0.81 −0.76 −1.54 5.25 −1.16 −1.17 0.28 5.25 −1.16 −1.17 0.28 −1.17 0.28 5.25 −1.16 −1.17 0.28 5.25 −1.16 5.25 −1.16 −1.17 0.28 −1.17 0.28 5.25 −1.16

25.31 3.097 25.48 450.247 2.67 553.354 16.713 -0.566 16.713 −0.566 2.675 553.355 3.173 −1.296 16.215 554.085 16.21 554.084 3.173 −1.296 16.431 442.931 10.628 108.052 17.048 106.613 10.011 444.371 10.628 108.053 16.431 442.931 10.011 444.371 17.048 106.612 4.654 107.189 −1.766 108.633 −1.766 108.629 4.654 107.193

020 011̅ ̅

01̅1

011̅

011

110

11̅0

1̅10

1̅1̅0

101

1̅01̅

a

The values in bold signify the most favorable docking position given by the lowest value of binding energy. All the energy values in the table are in kcal/mol. 1424



Article

Table 8. Native and Modified Slice and Attachment Energies of α-Glycine in the Presence of L-Alaninea slice

pos

02̅0

1 2 1 2 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

020 01̅1̅

01̅1

011̅

011

110

110̅

1̅10

1̅1̅0

101

10̅ 1̅

a

Esl

Eatt

Hhh

Esl′

ΔEatt ′′

−78.5 −78.5 −78.5 −78.5 −36.264 −36.264 −35.823 −35.823 −35.823 −35.823 −36.265 −36.265 −35.823 −35.823 −36.265 −36.265 −36.264 −36.264 −35.823 −35.823 −50.611 −50.611 −50.728 −50.728 −50.728 −50.728 −50.611 −50.611 −50.728 −50.728 −50.611 −50.611 −50.611 −50.611 −50.728 −50.728 −17.784 −17.784 −17.784 −17.784 −17.784 −17.784 −17.784 −17.784

−6.865 −6.865 −6.865 −6.865 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −24.725 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −17.415 −33.855 −33.855 −33.855 −33.855 −33.855 −33.855 −33.855 −33.855

−85.365 −85.365 −85.365 −85.365 −60.989 −60.989 −60.548 −60.548 −60.548 −60.548 −60.99 −60.99 −60.548 −60.548 −60.99 −60.99 −60.989 −60.989 −60.548 −60.548 −68.026 −68.026 −68.143 −68.143 −68.143 −68.143 −68.026 −68.026 −68.143 −68.143 −68.026 −68.026 −68.026 −68.026 −68.143 −68.143 −51.639 −51.639 −51.639 −51.639 −51.639 −51.639 −51.639 −51.639

−52.85 −75.313 −52.85 −75.313 −32.834 518.63 −31.89 −35.579 −31.89 −35.579 −32.83 518.63 −31.89 −35.579 −32.83 518.63 −32.834 518.63 −31.89 −35.579 −39.43 393.48 −38.93 57.044 −38.93 57.045 −39.43 393.48 −38.93 57.045 −39.43 393.48 −39.43 393.48 −38.93 57.044 −18.38 90.565 −18.38 90.569 −18.38 90.565 −18.38 90.569

−7.76 −6.892 −7.05 549.7 −25.3 −26.88 −10.14 −25.677 −10.14 −25.674 −25.31 −26.88 25.64 528.75 −9.53 −24.23 −9.53 −24.24 −25.64 528.76 −9.655 91.625 −18.47 428.92 −10.26 −18.465 −18.593 −17.6 −18.48 428.92 −9.65 91.63 −18.587 −17.603 −10.26 −18.46 −13.61 −34.49 −35.67 411.43 −35.71 411.43 −13.61 −34.494

Hhp −60.61 −82.205 −59.9 474.387 −58.134 491.75 −42.03 −61.256 −42.03 −61.253 −58.14 491.75 −6.25 493.171 −42.36 494.4 −42.364 494.39 −57.53 493.181 −49.085 485.105 −57.4 485.964 −49.19 38.58 −58.023 375.88 −57.41 485.965 −49.08 485.11 −58.017 375.877 −49.19 38.584 −31.99 56.075 −54.05 501.999 −54.09 501.995 −31.99 56.075

All the energy values in the table are in kcal/mol.

extent. It is a high-fidelity mechanistic model based on surface events. Figure 14 shows the dependence on the growth rate reduction of the (020̅ ) face of α-glycine as predicted by the model in the presence of different compositions of L-alanine. The additive composition on this (020̅ ) face at various weight fractions of the additive in the solution is given in Table 6. We see that as the composition of the additive increases, the growth rate of that face decreases to a greater extent, which is in line with expectations. Also, at low supersaturation, there is a significant dead zone (region in which the face does not grow at all). This dead zone is released at higher supersaturation. The width of the dead zone increases with an increase in additive concentration.

As discussed earlier, our experiments show a dispersion of crystal shapes. We see from Figure 14 that fluctuations in local supersaturation and in the local additive composition that a single crystal sees will lead to distributions in crystal shape. Also, as pointed out previously, we use a simple apparatus for our experiments in which it is very difficult to maintain constant supersaturation. Our crystallizer is a quiescent one, so it is not possible to prevent local variations in concentration. While performing experiments, we start with a very small amount of seed crystals (obtained from crystallizing the host in the absence of additive). However, we cannot prevent nucleation of new crystals in the presence of the additive when starting with these seed host crystals. Our model only predicts the change in morphology caused when starting with a native crystal in 1425



Article

Table 9. Native and Modified Slice and Attachment Energies of Paracetamol in the Presence of PAAa slice

pos

(110)

1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 1 2

(110̅ )

(1̅10)

(1̅1̅0)

(011)

(011̅)

(01̅1)

(01̅1)̅

(200) (2̅00) a

Esl −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −19.07 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −13.93 −11.42 −11.42 −11.42 −11.42

Eatt −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −6.56 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −9.135 −10.39 −10.39 −10.39 −10.39

Hhh −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −25.63 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −23.065 −21.81 −21.81 −21.81 −21.81

Esl′ ≈106 −19.2 ≈106 −19.3 −19.3 ≈106 −19.2 ≈106 −19.2 ≈106 −19.3 ≈106 ≈106 −19.3 ≈106 −19.2 ≈106 −14.15 ≈106 −14.15 −14.15 ≈106 −14.15 ≈106 −14.15 ≈106 −14.15 ≈106 ≈106 −14.15 ≈106 −14.15 −12.06 −12.06 −12.06 −12.06

Eatt ′′

Hhp

Eatt ′

−6.74 ≈106 −6.71 −6.43 −6.43 −6.71 ≈106 −6.74 ≈106 −6.74 −6.43 −6.71 −6.71 −6.43 −6.74 ≈106 −9.38 ≈106 ≈106 −9.26 ≈106 −9.38 −9.26 ≈106 −9.26 ≈106 ≈106 −9.38 ≈106 −9.26 −9.38 ≈106 −10.34 ≈106 ≈106 −10.34

≈106 ≈106 ≈106 −25.73 −25.73 ≈106 ≈106 ≈106 ≈106 ≈106 −25.73 ≈106 ≈106 −25.73 ≈106 ≈106 ≈106 ≈106 ≈106 −23.41 ≈106 ≈106 −23.41 ≈106 −23.41 ≈106 ≈106 ≈106 ≈106 −23.41 ≈106 ≈106 −22.4 ≈106 ≈106 −22.4

−6.52 −6.65 −6.69 −6.54 −6.48 −6.63 −6.63 −6.48 −6.65 −6.52 −6.55 −6.69 −6.62 −6.49 −6.49 −6.62 −9.08 −9.04 −9.08 −9.06 −9.03 −9.21 −9.2 −9.08 −9.2 −9.08 −9.03 −9.21 −9.08 −9.06 −9.08 −9.04 −10.39 −10.63 −10.63 −10.39

Δb ≈106 −0.22 ≈106 −0.21 −0.15 ≈106 −0.2 ≈106 −0.22 ≈106 −0.22 ≈106 ≈106 −0.16 ≈106 −0.19 ≈106 −0.125 ≈106 −0.145 −0.115 ≈106 −0.285 ≈106 −0.285 ≈106 −0.115 ≈106 ≈106 −0.145 ≈106 −0.125 −0.64 −0.88 −0.88 −0.64

All the energy values in the table are in kcal/mol. The values in bold signify the docking position in which the additive has the maximum effect.

a solution environment containing the additive; it does not describe the effect on nucleation. Hence, the experiments do not exactly mimic those conditions as there are seed as well as nucleated crystals in close proximity. We believe our model can be systematically improved to provide even better predictions of crystal morphology. The target points for increased fidelity are (1) Additive Uptake on Crystal Faces: Our current mean field model is adequate for predicting accurately which faces will recognize the additive and which will not. Our predictions are all made at low supersaturation, where the crystal surfaces are known to grow via spirals. (At higher supersaturation, it is common knowledge that the growth is mainly via 2-D clusters.) However, the mean field theory (eq 5) used for estimating the concentration of additive on a crystal surface has many assumptions which one needs to be careful about. For example, the model assumes that the molar volume of the host and the additive molecule are nearly equal, whereas in our case study, the molar volume of the PAA molecule is 1.26 times that of paracetamol. Also, the mean field equation assumes bulk thermodynamic quantities such as enthalpy of vaporization and enthalpy of solvation to represent the imposter interaction

energies, which may not always reflect the true chemical potential at the terrace-solvent interface. This is significant because the mean field equation is very sensitive to these interaction energy values / vv(l) and / vp(l) (because they appear inside the exponential term of the equation). It may be better to develop a more detailed method to estimate these interaction parameters rather than employing a first-pass method of using bulk thermodynamic quantities. Molecular simulations to estimate these interaction energy parameters might improve the model. We would further like to point out that though the dependence of critical length on S enables morphology prediction at different supersaturations, the mean field equation by itself is an equilibrium calculation (equilibrium distribution of imposter between a given crystal surface and solution) and additive composition on a crystal face depends only on the composition of the additive in solution and not on the composition of the host molecules. Chow et al.42 report in their experiments that the total additive uptake by the a paracetamol crystal changes significantly with an increase in supersaturation (solute concentration), shown in Figure 1c of their paper. Hence, the mean field theory is a simplification of some sorts and a nonequilibrium calculation or simulation that also accounts for a 1426



Article

images. We are also thankful to Rahul Sangodkar, UCSB, for his help in obtaining the SEM images.

positive constant chemical potential difference for the host molecules is a potential scope of improvement. The essential calculation would be analogous to the semigrand canonical ensemble described in Knott et al.47 and in Kofke and Glandt.48 Another point worth mentioning is that while calculating supersaturation, it is assumed that the additive does not greatly affect the solubility of the host. This assumption may not hold at higher additive concentrations where impurities are known to change the solubilities of the host molecules.49 (2) Solvent Effect for Modifying Kink Energies and Scaling of solid−solid interactions: The {201̅} faces of paracetamol disappear during the shape evolution in the presence of additive leading to bars (Figures 11 and 12). The experimental figures (Figure 13a−c) still show the small facets. A decrease of 10− 12% in the growth rate of {201̅} faces in the native morphology model will prevent these faces from disappearing from our predicted shapes. This is within the expected error when using a very simple model for solvent effect. An improvement in the solvent model by using either the solvent accessible surface area of the growth units at each kink site or a molecular scale estimate of the solvent effect rather than using the current method based on macroscopic solubility parameters will lead to a systematic improvement in morphology predictions. At present, it is not clear how to determine apriori whether a given edge will be inhibited by a reduction in step velocity (by either a kink density reduction or pinning of steps by immobile impurities) or by the reduction in spiral rotation time. It can be imagined that some faces exhibit step pinning whereas other faces exhibit spiral pinning on the same crystal.

■

(1) Variankaval, N.; Cote, A. S.; Doherty, M. F. AIChE J. 2008, 54, 1682−1688. (2) Yang, H. G.; Sun, C. H.; Qiao, S. Z.; Zou, J.; Liu, G.; Smith, S. C.; Cheng, H. M.; Lu, G. Q. Nature 2008, 453, 638−641. (3) Davis, K. J.; Dove, P. M.; De Yoreo, J. J. Science 2000, 290, 1134− 1137. (4) Coveney, P. V.; Davey, R.; Griffin, J. L. W.; He, Y.; Hamlin, J. D.; Stackhouse, S.; Whiting, A. J. Am. Chem. Soc. 2000, 122, 11557− 11558. (5) Weissbuch, I.; Addadi, L.; Berkovitch-Yellin, Z.; Gati, E.; Weinstein, S.; Lahav, M.; Leiserowitz, L. J. Am. Chem. Soc. 1983, 105, 6615−6621. (6) Weissbuch, I.; Addadi, L.; Lahav, M.; Leiserowitz, L. Science 1991, 253, 637−645. (7) Rimer, J. D.; An, Z.; Zhu, Z.; Lee, M. H.; Goldfarb, D. S.; Wesson, J. A.; Ward, M. D. Science 2010, 330, 337−341. (8) Clydesdale, G.; Roberts, K. J.; Lewitas, K.; Docherty, R. J. Cryst. Growth 1994, 141, 443−450. (9) Berkovitch-Yellin, Z. J. Am. Chem. Soc. 1985, 107, 8239−8253. (10) Chernov, A. Sov. Phys. Usp. 1961, 4, 117−145. (11) Cabrera, N.; Coleman, R. In The Art and Science of Growing Crystals; Gilman, J., Ed.; John Wiley & Sons, Inc.: New York, 1963. (12) Cabrera, N.; Vermilyea, D. In Growth and Perfection of Crystals; Doremus, R., Roberts, B., Turnbull, D., Eds.; John Wiley & Sons, Inc.: New York, 1958. (13) Sizemore, J. P.; Doherty, M. F. Cryst. Growth Des. 2009, 9, 2637−2645. (14) Ristic, R. I.; DeYoreo, J. J.; Chew, C. M. Cryst. Growth Des. 2008, 8, 1119−1122. (15) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. A 2009, 465, 1145− 1171. (16) Kaischew, R.; Budevski, E. Contemp. Phys. 1967, 8, 489−516. (17) Voronkov, V. V. Sov. Phys. Cryst. 1973, 18, 19−23. (18) Sizemore, J. P. Morphology Change of Molecular Organic Crystals via Molecular Imposters. Ph.D. Thesis, University of California, Santa Barbara, 2008. (19) Clydesdale, G.; Docherty, R.; Roberts, K. J. Comput. Phys. Commun. 1991, 64, 311−328. (20) Case, D. A. et al. AMBER 10; University of California: San Francisco, 2008. (21) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179−5197. (22) Frisch, M. J. et al. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2003. (23) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (24) McLachlan, A. D. Acta Crystallogr., Sect. A 1972, 28, 656−657. (25) Kabsch, W. Acta Crystallogr., Sect. A 1976, 32, 922−923. (26) Lesk, A. M. Acta Crystallogr., Sect. A 1986, 42, 110−113. (27) Thomas, T. N.; Land, T. A.; Martin, T.; Casey, W. H.; DeYoreo, J. J. J. Cryst. Growth 2004, 260, 566−579. (28) Lovette, M. A.; Doherty, M. F. J. Cryst. Growth 2011, 327, 117− 126. (29) Burton, W. K.; Cabrera, N.; Frank, F. C. Phil. Trans. R. Soc. A 1951, 243, 299−358. (30) Myerson, A. S.; Lo, P. Y. J. Cryst. Growth 1990, 99, 1048−1052. (31) Ward, M. D. Chem. Rev. 2001, 101, 1697−1725. (32) Bisker-Leib, V.; Doherty, M. Cryst. Growth Des. 2003, 3, 221− 237. (33) Boek, E. S.; Feil, D.; Briels, W. J.; Bennema, P. J. Cryst. Growth 1991, 114, 389−410. (34) Winn, D.; Doherty, M. F. AIChE J. 1998, 44, 2501−2514. (35) Winn, D.; Doherty, M. F. AIChE J. 2000, 46, 1348−1367. (36) Barton, A. F. M. Chem. Rev. 1975, 75, 731−753.

7. CONCLUSIONS We have developed a mechanistic model that can predict the crystal morphology of any molecular organic crystal, in the presence of a structurally similar additive. The model explains well the morphology change of α-glycine in the presence of L-alanine, as well as paracetamol in the presence of PAA. This method is many orders of magnitude faster than carrying out all-atom simulations for predicting crystal growth rates and the effect of additives on crystal surfaces. This model and the generic understanding it provides of the mechanistic inhibition of an additive on crystal growth is expected to be very useful for screening or designing tailor-made additives to engineer crystal morphology. The model, in its current form, is not very well suited for macromolecular additives such as polymers, but the underlying mechanistic principles of disruption of either the step velocity or the critical length are still expected to be the same for any kind of additive. Detailed in situ atomic force microscopy experiments might reveal which step edges actually exhibit the spiral pinning mechanism as opposed to the step velocity reduction mechanism. This is a subject of ongoing research.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*Phone: 805-893-5309. Fax: 805-893-4731. E-mail: mfd@ engineering.ucsb.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful for the financial support provided by Abbott Laboratories and Eli Lilly. We are thankful to Daniel Markarian for his assistance with the experiments and capturing optical 1427



Article

(37) Gaffney, J. S.; Pierce, R. C.; Friedman, L. J. Am. Chem. Soc. 1977, 99, 4293−4298. (38) Kuvadia, Z. B.; Doherty, M. F. Cryst. Growth Des. 2011, 11, 2780−2802. (39) Veesler, S.; Lafferrere, L.; Garcia, E.; Hoff, C. Org. Process Res. Dev. 2003, 7, 983−989. (40) Lovette, M. A.; Muratore, M.; Doherty, M. F. AIChE J. 2012, 58, 1465−1474. (41) Ristic, R.; Finnie, S.; Sheen, D.; Sherwood, J. J. Phys. Chem. B 2001, 105, 9057−9066. (42) Chow, A.-L.; Chow, P.; Zhongshan, W.; Grant, D. Int. J. Pharm. 1985, 24, 239−258. (43) Kuvadia, Z. B. Crystal Morphology Prediction and Modification of Organic Molecular Crystals. Ph.D. thesis, University of California, Santa Barbara, 2012. (44) Granberg, R. A.; Rasmuson, k. C. J. Chem. Eng. Data 1999, 44, 1391−1395. (45) Hendriksen, B. A.; Grant, D. J.; Meenan, P.; Green, D. A. J. Cryst. Growth 1998, 183, 629−640. (46) Zhang, Y.; Sizemore, J. P.; Doherty, M. F. AIChE J. 2006, 52, 1906−1915. (47) Knott, B. C.; Doherty, M. F.; Peters, B. J. Chem. Phys. 2011, 134, 154501. (48) Kofke, D. A.; Glandt, E. D. Mol. Phys. 1988, 64, 1105−1131. (49) Sangwal, K. Additives and Crystallization Processes: From Fundamentals to Applications; John Wiley & Sons, Ltd: New York, 2007.

1428


Effect of Structurally Similar Additives on Crystal Habit of Organic

Recommend Documents