Spiral Growth Model for Faceted Crystals of Non ... - ACS Publications

ARTICLE pubs.acs.org/crystal

Spiral Growth Model for Faceted Crystals of Non-Centrosymmetric Organic Molecules Grown from Solution Zubin B. Kuvadia and Michael F. Doherty* Department of Chemical Engineering, University of California Santa Barbara, Santa Barbara, California 93106-5080, United States

bS Supporting Information ABSTRACT: The classical BurtonCabreraFrank (BCF) spiral growth model fails to work satisfactorily for many noncentrosymmetric organic molecules such as active pharmaceutical ingredients (APIs), nonlinear optical compounds, etc., due to the inherent assumption of Kossel crystal structure in the solid-state. We develop a more general mechanistic spiral growth model that enables morphology prediction for all kinds of organic molecules. We develop generalized expressions for kink density, kink incorporation rate, and step velocities for such molecules. Stable and unstable edges arising due to the complex bonding pattern are discussed and treated. We demonstrate the applicability of the model by correctly predicting the experimental morphologies of paracetamol and lovastatin grown from solution.

1. INTRODUCTION AND MOTIVATION Crystallization from solution is a vital unit operation for the industrial production of various chemicals in solid form, for example, pharmaceutical products, catalysts, specialty and bulk chemicals, and chemical intermediates. This process defines many fundamental properties of the resultant materials including chemical purity and composition, polymorphic state, crystal size, shape distributions, etc. Crystal morphology is known to significantly influence the end-use efficacy of solid products (e.g., bioavailability for pharmaceutical compounds,1 reactivity for catalysts2), as well as the downstream performance of the entire process (e.g., by affecting filtering and drying times1). The desired shape and size distributions of a crystalline product depend on the specific application; hence, a general methodology for the prediction and improvement of crystal morphology would serve as a major breakthrough for product and process design. Our main focus is on molecular organic systems (of pharmaceutical interest or other applications such as nonlinear optics, organic semiconductors, etc), with a view to the development and improvement of fast, accurate models based on firstprinciples that predict the crystal morphology. 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 a connection between the surface and the bulk solid-state chemistry. It is well established by detailed experiments36 that, at low supersaturation, crystals of organic molecules grow via spirals emanating from screw dislocations. The spiral growth model of Burton, Cabrera, and Frank,7 widely recognized as the BCF model, was developed six decades ago, and several modifications and extensions exist, including Chernov’s work8 and the spiral growth model of Snyder and Doherty.9 In the work by Snyder r 2011 American Chemical Society

and Doherty, the systematic implementation of the BCF model was successfully demonstrated for organic molecules or crystal growth units having a center of symmetry (centrosymmetric). It will be shown in Section 4 of this article that growth units with a center of symmetry behave in an isotropic fashion. However, the vast majority of molecules of realistic complexity do not have a center of symmetry, and they warrant a much more rigorous treatment of the spiral growth model to give accurate crystal shape prediction. In the past two decades, there have been several important developments in the field of non-Kossel crystal growth. Zhang and Nancollas worked on the step movement of AB ionic crystals having a NaCl lattice.10 Chernov et al.11,12 derived expressions for the step velocity of non-Kossel crystals, mainly addressing the surface physics and validating them for AB type ionic systems such as CaCO3. Cuppen et al.13,14 treated molecular Non-Kossel crystals of two asymmetric growth units in the crystal structure using Monte Carlo simulations and discussed kink incorporation and step growth for a specific crystal bonding structure derived from naphthalene (which itself is centrosymmetric). In this work, we describe a general mechanistic model to predict the crystal morphology of any organic molecular crystal (centrosymmetric as well as non-centrosymmetric) belonging to any space group and having any number of asymmetric units within the unit cell. The only inputs required are the crystal structure, a force-field model capable of accurately representing the solid-state chemistry, and solvent information. We address the problem from a mechanistic perspective (in contrast to an Received: November 23, 2010 Revised: March 18, 2011 Published: March 24, 2011 2780

dx.doi.org/10.1021/cg101560u | Cryst. Growth Des. 2011, 11, 2780–2802

Crystal Growth & Design

ARTICLE

unconstrained molecular dynamics/quantum mechanical viewpoint) and outline the method to solve the problem from start to finish. This paper is organized in the following way: first we give a brief background of the traditional spiral growth model for predicting crystal shape, then we define centrosymmetric and non-centrosymmetric growth units and their relevance in the context of crystal growth and highlight the differences in crystal growth between them. We then present a model to calculate step velocity and hence predict morphology for crystals of noncentrosymmetric growth units using paracetamol as a central example. Finally, we apply the theory developed to two compounds of pharmaceutical interest, paracetamol and lovastatin, and compare the predictions with experimental data from the literature and from our laboratory.

2. THEORETICAL BACKGROUND OF THE SPIRAL GROWTH MODEL The growth rate of a crystal face that is growing by the spiral mechanism can be expressed as h τ

ð1Þ

Figure 1. Examples of centrosymmetric growth units. (a) Naphthalene, (b) adipic acid.

where τ is the characteristic rotation time of the spiral, which corresponds to the time between consecutive passes of steps across a point on the face, and h is the step height.9 Snyder and Doherty9 give an expression for the rotation time of a convex Nsided polygonal spiral as

dynamic shape using the faceted evolution model developed by Doherty et al.9,23

G¼

τ¼

N

l

∑ c, i þ 1 i¼1

sinðRi, i þ 1 Þ vi

ð2Þ

where Ri,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 Budevski15 and later by Voronkov16 v ¼ 0 l e lc v ¼ v¥ l > lc

ð3Þ

whereby a spiral edge starts growing with a constant step velocity only after reaching a critical length. The concept of a critical length has also been experimentally validated for calcite17 and for lysozyme.11 Given the crystallography of the system under consideration and using modern force fields like Assisted Model Building with Energy Refinement (AMBER)18,19 to calculate intermolecular interactions, it is possible to determine the periodic bond chains (PBCs) using the rules proposed by Hartman and Perdok.20 These PBCs comprise chains of strong interactions and are taken to be the sides of the spiral and thus the direction of the steps flowing across the surface. We will explain in detail the methods to calculate step velocity for crystal systems having centrosymmetric growth units in Section 4 and noncentrosymmetric growth units in Section 5.2. This facilitates the estimation of relative growth rates of the relevant faces of the crystal and the construction of either the steady-state shape using the FrankChernov condition21,22 G1 G2 Gi ¼ ¼ 333 ¼ H1 H2 Hi

ð4Þ

where Gi is the perpendicular growth rate of face i and Hi is the perpendicular distance of that face from the crystal center, or the

3. KOSSEL, NON-KOSSEL, CENTROSYMMETRIC, AND NON-CENTROSYMMETRIC GROWTH UNITS It is worthwhile reviewing what exactly the terms Kossel and non-Kossel crystal systems mean and also to remind the reader about the concepts of centrosymmetric and non-centrosymmetric growth units in the context of crystal growth and morphology prediction. (Symmetry of the center of a molecule has been used commonly in solid state chemistry but not so extensively in the specific field of crystal growth.) A Kossel crystal is a simple cubic lattice structure with equal bond strength in all (six) directions.24 The only known substance that satisfies these requirements and hence qualifies as a Kossel crystal is one of the polymorphs of the metal polonium.25 All other elements and molecules are Non-Kossel and can be further classified as centrosymmetric or non-centrosymmetric based on the arrangement of individual atoms within the molecule. Centrosymmetric molecules are molecules in which the center of mass is also an inversion center. In such a molecule, for every point (x, y, z), there is an indistinguishable point (x, y, z) (e.g., the number 8 is an example of a centrosymmetric object). Naphthalene and adipic acid are examples of centrosymmetric molecules; see Figure 1. Non-centrosymmetric growth units do not have a center of inversion, for example, acetaminophen (also commonly referred as paracetamol); see Figure 2. In crystal structures where the growth unit is more than a single molecule (a dimer, tetramer etc.), the center of symmetry of the crystallographic growth unit is calculated with respect to the entire growth unit, rather than the single molecule. As seen in Figure 3a,b, a single molecule of ibuprofen does not have a center of symmetry, but the dimer is centrosymmetric. The crystal growth of ibuprofen is known to occur by addition of dimers to the lattice; hence ibuprofen crystals are considered to grow by addition of centrosymmetric growth units. The hydrogen bonded 2781

dx.doi.org/10.1021/cg101560u |Cryst. Growth Des. 2011, 11, 2780–2802

Crystal Growth & Design dimers of ibuprofen molecules with the centers of mass of the dimers indicated by a dot are shown in Figure 3c. It is important to note the distinction between the symmetry of the crystallographic growth unit and the symmetry of the crystal lattice as a whole (centrosymmetric and non-centrosymmetric space groups). A centrosymmetric crystal system refers to a space group which contains an inversion center as one of its symmetry elements. In such a space group, for every point (x, y, z) in the unit cell there is an indistinguishable point (x, y, z). P21/c, P2/c, and Pbca are examples of centrosymmetric space groups, while P212121 and P21 are examples of non-centrosymmetric space groups. Centrosymmetric growth units can occur in centrosymmetric as well as non-centrosymmetric space groups and vice versa.26 In this study, we focus on centrosymmetric and non-centrosymmetric molecules/growth units in all possible space groups, as the symmetry of the crystallographic growth unit dictates whether the system acts in a Kossel-like or a Non-Kossel fashion. Crystals of centrosymmetric growth units, though non-Kossel by definition, do not deviate significantly from Kossel crystals and can be considered to behave in a Kossel-like fashion (represented by green color in Table 1) irrespective of whether they are forming in centrosymmetric or non-centrosymmetric space groups. The behavior of centrosymmetric growth units is discussed in the following section. Non-centrosymmetric growth units lead to complex periodic bonding networks and behave in a

ARTICLE

Non-Kossel manner (represented by the red color in Table 1). Again, this behavior is independent of whether the space group of the crystal is centrosymmetric or non-centrosymmetric. Our analysis and growth rate model applies to all types of growth units (centrosymmetric and non-centrosymmetric) crystallizing in all types of space groups. Therefore, it has wide applicability to crystalline products of realistic complexity, for example, APIs, dyes, nonlinear optical materials, protein crystals, etc.

4. DIFFERENCES BETWEEN CENTROSYMMETRIC AND NON-CENTROSYMMETRIC GROWTH UNITS: CHALLENGES FOR MODELING NON-CENTROSYMMETRIC GROWTH UNITS The main difference between crystals of centrosymmetric and non-centrosymmmetric growth units is the isotropic driving force experienced by centrosymmetric growth units to incorporate into the solid state structure. Consider the bonding structure for the crystal system of naphthalene (a centrosymmetric molecule) with Figure 4 showing three different faces. Naphthalene molecules within the unit cell have two different symmetry-related orientations possible, shown by green and black color dots that represent the centers of mass of these growth units. There are six different bond directions and hence six bond chains possible for the entire crystal system. (All the periodic bond chains in this work were drawn using computerized programs for visualizing bond chains developed by Lovette and Doherty.27) Table 1. Examples of Centrosymmetric and Non-Centrosymmetric Growth Units and Space Groups

Figure 2. Example of non-centrosymmetric growth unit, acetaminophen.

Figure 3. (a) Ibuprofen monomer, non-centrosymmetric, (b) ibuprofen dimer, centrosymmetric, (c) crystal packing structure of ibuprofen in the form of dimers with their centers of mass indicated by a dot. 2782



ARTICLE

Figure 4. Bonding structure for different faces of naphthalene, (a) {001} face, (b) {201} face, (c) {110} face. (A bond chain with an asterisk identifier is related to the bond chain of the same letter by symmetry.)

Figure 5. Centrosymmetric growth units behave in a Kossel-like fashion; that is, all kink sites are equivalent as the energy exposed at all kink sites is the same. (a) Bonding structure of {110} face of naphthalene with three periodic bond chains showing kink sites and bonds exposed at kink sites. (b) Bonding structure of {001} face of naphthalene with three periodic bond chains showing a kink site.

The selected {110} face has three independent periodic bond chains (h*, j*, and k; see Figure 5a). Since the growth unit is centrosymmetric, the bonds are equal in strength along a single direction. If we move along the bond chain direction h* in Figure 5a and consider kink sites (represented by gray dots), the bonds exposed at any kink site along this edge are the same (bonds φh*, φj*, and φk in the slice and the remaining three bonds hidden from view that are pointing upward from the layer below, φh, φj, and φi). The energy exposed at each kink site is half the lattice energy (Elatt = 2*(φh þ φh* þ φi þ φj þ φj* þ φk)). Hence, every kink site is a half-crystal position.7,24 Similarly, a kink site along any other periodic bond chain direction on this face (j* or k) also exposes the same six bonds. Hence, all kink sites are equivalent on the face. Hence, the driving force for the incorporation of solute molecules into the kink sites on this face is the same irrespective of the edge direction. This concept extends to all other faces for this centrosymmetric bonding structure as

shown in Figure 5b for the {001} face; the bonds in the slice are different (φh, φh*, φi), however, the net bonds exposed at the kink site are still the same six bonds. Irrespective of the edge or the face, there is just one type of kink site for the whole crystal lattice, which exposes the same bonds at the kink site (half the total number of bonds and energy exposed that is half the lattice energy). This means that the driving force for crystal growth is isotropic, implying it will have an equal effect on all orientations and hence the isotropic terms are not retained in the equations for calculating the relative growth rates for faces. The velocity of a step/edge is dependent only on the number density of kink sites on an edge and the geometry of the crystal. Hence, the crystal of a centrosymmetric growth unit, although being Non-Kossel, behaves in a Kossel-like fashion. An expression for the step velocity of the ith edge where kink sites along step edges are the primary sites for the incorporation of centrosymmetric growth units into the crystal and disincorporation 2783



ARTICLE

Figure 6. (a) Bonding structure of {110} face of naphthalene (centrosymmetric growth units), (b) bonding structure of {011} face of paracetamol (non-centrosymmetric growth units).

"

from the crystal, and where kink integration is the rate-limiting event is given by ae ΔU vi ¼ ap, i ν0 exp Vm ðC Csat Þ kT δ i

ð6Þ

where ap,i is the distance an edge i propagates with the addition of a single row of growth units (Figure 6a) and v0 is the frequency of attachment and detachment attempts.25 The average spacing between kinks is δ and the molecular spacing along the step is ae; thus, the quantity ae/δ is the probability of finding a kink at any location along the step edge; in other words, a measure of the density of kink sites (Fi) on the particular edge. C and Csat are the actual and the equilibrium concentrations of the solute in solution near the edge that is growing (with S = C/Csat, the supersaturation ratio) and Vm is the molar volume of the solute. Finally, ΔU is the energy barrier for the incorporation of molecules into a kink site. For solution crystallization, this quantity can be estimated as the enthalpy of desolvation. The frequency (v0), the difference between the actual and saturation concentration (C Csat), and the exponential (exp((ΔU/kT)) terms are all taken to be isotropic for centrosymmetric growth units as they exhibit Kossel-like behavior, and all these terms drop out in the calculation of relative growth rates for these crystals and therefore have no effect on crystal shape. The final form of the step velocity is then given by vi ap, i Fi

!#1 ð8Þ

ð5Þ

which becomes ae ΔU vi ¼ ap, i ν0 exp Vm Csat ðS 1Þ kT δ i

φkink vi ap, i 1 þ 0:5 exp RT

ð7Þ

We assume that the rearrangement of kinks occurs on a time scale that is faster than that of kink integration;28 therefore, the kinks are assumed to be in their most probable distribution and are Boltzmann distributed. This leads to

The term in brackets is the expression for kink density or the probability of finding a kink site in a step, and φkink is the work or the energy change required to create one kink site along an edge by rearrangement of the structure.7,9,29 Having established that crystals of centrosymmetric molecules behave in a Kossel-like fashion, we proceed to identify the specific challenges for noncentrosymmetric growth units. The major ones are listed as follows: (1) Complex PBCs: The periodic bond chains are much more complex in nature and additional rules must be imposed to find the periodic bond chains that constitute the most stable edges that form the sides of the spirals. Figure 6a shows the bonding structure for the {110} face of naphthalene (centrosymmetric). We select one chain [110] (red color) and moving along this chain, we see that all bonds that comprise this bond chain are equal in strength. Also, the chain divides the remaining bonds into two halves. The number of bonds and the strength of the bonds are equal in each half. All the growth units that make up this chain have the same set of bonds (2 red, 2 blue, and 2 yellow that fall within the slice). Therefore, even though the molecules attach to the crystal phase with different symmetry operators within the lattice, since the set of bonds that they satisfy when they incorporate is the same, the growth unit is essentially of just one kind. These attributes of the bonding structure for naphthalene are in contrast to those for paracetamol which is a non-centrosymmetric molecule for which the growth unit is the monomer. Considering the bonding structure for one of the faces of paracetamol (Figure 6b), the bonds comprising a single bond chain are not necessarily of the same strength. Also, the broken bond energy of each half of the slice divided by a bond chain are not necessarily equal. The adjacent molecules (with different symmetry operators) on 2784



ARTICLE

a chain may not have the same set of bonds due to the complex bonding, and hence they behave as different growth units forming different types of kink sites. The maximum number of different growth units for a particular crystal system is the same as the total number of asymmetric units in the unit cell. (2) Kink Rates: As discussed above, centrosymmetric growth units experience an isotropic driving force. Non-centrosymmetric growth units have more than one type of kink site and each edge on each face differs in terms of how many different types of kink sites it can have. Hence, the driving force for the incorporation of solute molecules into the kink sites is no longer isotropic. The velocity of the step needs to have a specific term which describes the net rate at which solute molecules incorporate within the different types of kink sites. The nature of this term will depend on the number and nature of the different types of kink sites on the particular edge. (3) Stable and Unstable Edges: Because of the complex bonding structure for non-centrosymmetric growth units and the fact that the slice may not be divided into two identical halves by a given PBC, different layers of that edge direction may be stable or unstable. This concept of unstable edges and how to treat them for the estimation of relative growth rate of faces is discussed in Section 5.3. Centrosymmetric growth units can only lead to stable edges due to the symmetric bonding structure.

crystal growth of complex noncentrosymmetric organic molecules. Paracetamol has been widely studied in the literature because of its significance in the pharmaceutical industry. There have been various attempts to study the shape evolution and develop models to predict its crystal shape. While experimental shape studies can be found in the literature (the most prominent among them being work by Ristic et al.30 and Shekunov and Grant31), the theoretical models to explain the particular shape behavior using the attachment energy model or kinetic Monte Carlo simulations have been, at best, only partially successful.3235 A complete understanding is still missing even after several years of focused attempts by several research groups all over the world. This arises mainly from the challenge in understanding complex bond chain networks with multiple types of kink sites applied to crystal growth and over reliance on the traditional Kossel model. Paracetamol crystallizes in three different polymorphic forms: a monoclinic form, which is shown in Figure 7; an orthorhombic form, and a trihydrate. Since the monoclinic form is the stable polymorph and pharmaceutically most relevant, we only discuss this form grown from aqueous solution. 5.1. Crystal Structure, Periodic Bond Chains, and Spiral Edges. Monoclinic paracetamol crystallizes in the space group P21/a. There are four molecules per unit cell related to each other by symmetry operations; see Figure 7. The cell parameters are a = 12.92 Å, b = 9.40 Å, c = 7.10 Å, and β = 115.9 as measured by Haisa et al.36 at room temperature. To predict the growth rate of different faces of a paracetamol crystal, the first step is to calculate the energy of interaction between the paracetamol molecules in the solid state. This was carried out using the program HABIT37 and applying the AMBER Force field18,19 as implemented by our group at UC Santa Barbara.38 The partial atomic charges on each of the atoms of the four molecules comprising the unit cell were calculated using GAUSSIAN 0339 with the restrained electrostatic potential (RESP) model.40 The different atomatom interactions get condensed into a single intermolecular interaction leading to a bond chain as shown in Figure 8. Table 2 shows the different bonds in the bonding structure for paracetamol calculated using the AMBER force field. The calculated value of the lattice energy of paracetamol is 28.27 kcal/mol, which compares remarkably well with the experimental value of 28.18 kcal/mol, measured by Li and Feng.41 To understand the bonding structure, two different views perpendicular to the c and b axes, respectively, are presented in Figure 9a,b and bonds of different strengths are shown using different colors. Also, the four different growth units related by symmetry operators that make up the unit cell are shown using different colors. Two of the molecules (labeled 1 and 4) are part of a two-dimensional (2-D) hydrogen bonded sheet perpendicular to the b-axis (see Figure 9b). The other two

5. SPIRAL GROWTH MODEL FOR CRYSTALS OF NONCENTROSYMMETRIC ORGANIC GROWTH UNITS Acetaminophen (paracetamol) is chosen to be the central example to demonstrate the resolution of the conundrum of

Figure 7. Paracetamol unit cell with four growth units (Z = 4) having different symmetry operators.

Table 2. Paracetamol Solid-State Interactions bond direction

bond distance (Å)

total bond energy (kcal/mol)

dispersive energy (kcal/mol)

coulombic energy (kcal/mol)

[100]

6.65

2.72

0.85

1.87

[010]

3.83

4.45

3.44

1.01 1.20

[010]

6.77

2.43

1.23

[102]

7.41

3.14

0.20

3.34

[112]

7.04

2.39

1.20

1.19

[011] [001]

5.90 7.08

1.06 1.38

1.01 0.41

0.50 0.97

[110]

8.31

0.50

0.53

0.03

2785



ARTICLE

Figure 8. (a) Various atomatom interactions between paracetamol molecules along the [100] direction. (b) The atomatom interactions represented by a single bond chain and molecules replaced by their centers of mass.

Figure 9. Paracetamol bonding structure; two different views (a) {001} face, (b) two parallel {020} planes. 2786



ARTICLE

Figure 10. Paracetamol bonding structure of {011} face; multiple intersecting periodic bond chains. (a) Bonding structure, (b) PBCs a a and b b forming the spiral edges, (c) PBCs a a and c c forming the spiral edges.

molecules (labeled 2 and 3) are part of a parallel sheet, translated by half a unit cell down the b-axis. Within each sheet there are two hydrogen bond chains (four H-bonds per molecule illustrated by red and yellow colored bond chains). The non-centrosymmetric nature of the bonding network is highlighted when one moves along the [010] chain in Figure 9a where we see bonds of unequal strength (shown in black and green colors) in a single direction. The periodic bond chains (PBCs) were constructed from the strong repeating intermolecular bonds. From the complex bond network, only the straight chains that can act as a ledge (step) for formation of kinks and incorporation of growth units in it were selected. On each face, we need to determine the edge energy and kink energy for each bond chain to be able to calculate the step velocity. Hartman and Perdok42 state some rules for selecting the PBCs in their paper,42 “Chains of strong bonds are sought. It must be possible to divide the structure of an F-face into a fabric of chains that do not have any bond in common. The chain should not entirely consist of periods of other chains.” Considering the {001} face in Figure 9a, it can be seen that the molecules fall in a straight chain in two different directions; the horizontal ([100] edge) and the vertical ([010] edge) directions. (No other chain can be constructed without violating the condition of not having any common bonds.) This means that a four-sided spiral will be formed on this face. Non-centrosymmetric molecules generally lead to complex bond chains; hence for some faces (like the {200}, {011}, and {110} faces of paracetamol), multiple intersecting chains can be formed for a single face from which we need to select the most stable PBCs that can form kink sites for the incorporation of solute molecules. These PBCs should form a closed polygon which then correspond to the edges of a faceted spiral and thus enables the face to grow layer-by-layer. For example, consider the {011} face. The bonding network is shown in Figure 10. A closed polygon can be formed from bond chains a and b, as shown in

Figure 10b or an alternate polygon comprising bond chains a and c can be formed as shown in Figure 10c. The bond chains b and c, however, are mutually exclusive since they share a common bond (green-colored bond); that is, these two chains cannot coexist to form one spiral polygon. Hence, we select the stronger bond chain (the bonds that comprise the chain should be of higher strength) as this will lead to a more stable edge and hence a more stable faceted spiral.43 For the {011} face of paracetamol, bond chain b is stronger along the edge (average energy per bond, 2.915 kcal/mol) than bond chain c (average energy per bond, 1.745 kcal/mol); hence it is selected. Using the calculations of the intermolecular bond strengths and the periodic bond chains, we identified the families of F-faces of paracetamol (the faces that have two or more periodic bond chains and can, therefore, form spiral growth rings and grow layer-by-layer) to be {001}, {201}, {110}, {011}, {200}, and {020}; see Figure 11. 5.2. Step Velocity for Non-Centrosymmetric Organic Growth Units. The simplifying assumption of isotropic behavior on all step fronts that made centrosymmetric growth units a special case is not applicable for non-centrosymmetric growth units. We need to derive a modified expression for step velocity that will be valid for all organic growth units and also simplify to the known expressions (eqs 5, 6, and 7) for centrosymmetric growth units. Since kink sites are the incorporation sites for the solute molecules, the step velocity is the product of the distance of propagation, density of the kink sites, and the rate of incorporation into the kink sites. Hence, vi ap , i u i Fi

ð9Þ

The kink rate term, ui, now accounts for the nonisotropic behavior of the non-centrosymmetric growth units (for centrosymmetric growth units u is independent of the specific edge and face and hence is the same for all step fronts on all faces). This 2787



ARTICLE

Figure 11. Bonding structure of the different F-faces of paracetamol (a) {001} face, (b) {011} face, (c) {200} face, (d) {201} face, (e) {110} face, (f) {020} face.

Figure 12. Types of kink sites on the [010] edge of the {001} face of paracetamol. (a) The [010] straight edge, (b) four different types of kink sites on the [010] edge based on the bonds exposed at the kink sites.

approach of using the kink rate to estimate the step velocity has been adopted and justified by Chernov et al.,11,44 Zhang and Nancollas,10 and Cuppen et al.13,14 to describe Non-Kossel crystal growth. These researchers have shown how to develop expressions for ui when the growth units number no more than three. However, paracetamol and lovastatin each have four growth units; RDX has eight,45 and larger numbers may be expected for more complex systems. In the following sections, we

will demonstrate how to calculate kink density and kink rates on selected edges of relevant faces of paracetamol and also explain how the concepts can be generalized for all edges and faces of any organic crystal belonging to any space group. 5.2.1. Determination of Kink Densities. Since spiral edges on the faces of a crystal of non-centrosymmetric growth units have more than one type of kink site, it is convenient to group and define the kink sites based on the specific bonds they expose. 2788



ARTICLE

Figure 13. Nomenclature of the types of kink sites on the [010] edge of the {001} face of paracetamol based on the terminating growth unit and orientation, (a) Af, (b) A), (c) Bf, and (d) B).

Consider the [010] edge of the {001} face of paracetamol. A pictorial representation of the edge and the different kink sites is provided in Figure 12. The selected straight edge [010] has two different growth units (shown by turquoise and black dots). Though the growth units are the same chemical entity and indistinguishable in solution, the symmetry and orientation with which they attach to the solid phase are different and this symmetry/orientation defines the growth unit. From Figure 12b, we see that four different types of kink sites, Af, A), Bf, and B), can exist on this edge. The kink sites are defined by the bonds that they expose (each of the kink sites expose different bonds); the way that these kink sites are named is as follows. If the segment of the leading edge (that breaks to form a kink) terminates with the first type of growth unit (turquoise colored dot), the kink site is labeled as A, and likewise if the leading edge terminates with the second type

of growth unit (black dot), the kink site is labeled as B type. Further, the kink site of type A itself can be of two different types depending on whether the segment of the edge terminates from the left or the right, as shown in Figure 13a,b. (Same is the case with kink site of type B-see Figure 13c,d. We name these kink sites of the type f and ) arbitrarily (choosing either orientation as f, the other one will be ).) We could have used “left” and “right” to name the kinks, but we prefer the f and ) notation as it is independent of the view. Hence, for two growth units on this edge, we end up with four different types of kink sites. Close to equilibrium (i.e., at low supersaturation), the kinks are assumed to form by thermal rearrangement of a straight step. We still follow the BCF approach of finding kink density.7 However, since we have multiple types of kinks and these kinks occur in different pairs with different energies of formation (energy required to form kinks from a straight step by the 2789



ARTICLE

Figure 14. Different rearrangements of the edge to form different types of kink sites on the [010] edge of the {001} face of paracetamol. (a) Rearrangement type I ((3,4) f (1,2) transition) forms 2 Af and 2 B) kink sites. (b) Rearrangement type II ((3,4) f (2,1) transition) forms 1 Af, 1 B), 1 Bf, and 1 A) kink sites. (c) Rearrangement type III ((4,3) f (1,2) transition) forms 1 Af, 1 B), 1 Bf, and 1 A) kink sites. (d) Rearrangement type IV((4,3) f (2,1) transition) forms 2 Bf and 2 A) kink sites.

Figure 15. Determination of kink rates. (a) Incorporation of a growth unit into a kink site and a hypothetical step with four different types of kink sites in series. (b) Steady-state 1D master equation to find probability of different types of kink sites.

process of rearrangement), we cannot use the classical Boltzmann formula for kink density as we did in eq 8. Instead, we count all the different microstates (ways of rearrangement from a

straight edge) that would give each of the kink site types and use the Boltzmann distribution to find each of their probabilities. The sum of all their probabilities provides the total number density of 2790



ARTICLE

Figure 16. Stable and unstable edges. The [001] edge of the {200} face of paracetamol; a) top view of the face, b) stable edge, c) unstable edge.

all the kink sites. Consider the Figure 14ad. These figures show the four different possible ways to rearrange the [010] edge and form kinks on it. Figure 14a outlines one of these possible rearrangements, labeled as I ((3,4) f (1,2) transition) which forms two Af kink sites and two B) kinks. The energy expense for this thermal rearrangement per kink formed is ε1 kcal/mol (ε1 = (4.45 þ 4.45)/4)). Likewise, the Figure 14bd depict the rearrangements II ((3,4) f (2,1) transition), III ((4,3) f (1,2) transition) and IV ((4,3) f (2,1) transition); the types of kinks formed and the energy expense per kink, ε2 = (4.45 þ 2.43)/4, ε4 = (4.45 þ 2.43)/4, and ε4 = (4.45 þ 2.43)/4, respectively. Keeping track of the number of kinks of each type (degeneracy) (for example, considering Af kinks; two are formed by rearrangement of type I, one by rearrangement of type II and one by type III) and knowing the energy expense for each rearrangement, we can determine the probability of formation for each type of kink site and add up all the individual probabilities to find the total kink density, F, on this edge. jε1 j jε2 j jε3 j þ exp þ exp =ð4Q Þ 2 exp kT kT kT jε4 j jε2 j jε3 j FBf ¼ 2 exp þ exp þ exp =ð4Q Þ kT kT kT jε4 j jε2 j jε3 j FA) ¼ 2 exp þ exp þ exp =ð4Q Þ kT kT kT jε1 j jε2 j jε3 j þ exp þ exp =ð4Q Þ FB) ¼ 2 exp kT kT kT FAf ¼

ð10Þ Therefore, F ¼ FAf þ FBf þ FA) þ FB)

ð11Þ

giving, 1 jε1 j jε2 j jε3 j jε4 j þ exp þ exp þ exp exp C B C B kT kT kT kT C F¼B @ jε1 j jε2 j jε3 j jε4 j A 1 þ exp þ exp þ exp þ exp kT kT kT kT 0

ð12Þ where Q is the partition function of the Boltzmann distribution, given by

jε1 j jε2 j jε3 j jε4 j þ exp þ exp þ exp Q ¼ 1 þ exp kT kT kT kT

The leading term in Q accounts for the probability that there is no kink at a site. The same sequence of steps can be repeated to find the kink density on any step of any face. 5.2.2. Determination of Kink Rates. Since crystal growth is a kinetic process, the net kink rate ui is the difference between two opposing fluxes; the attachment rate of solute molecules into the kink sites and the detachment rate of solute molecules from the kink sites. The attachment rate of the solute into the kink, R, is considered independent of the exact position of attachment; it is only dependent on the driving force (S is the supersaturation ratio of C/Csat and v0 is the frequency of attachment and detachment attempts), whereas the detachment rate from a particular kink site, vk,i, depends on the type of kink site and the specific bonds being broken.10,46 Elatt R ¼ ν0 S exp 2kT 0 1 ðφk, i Þj ð13Þ j B C νk, i ¼ ν0 exp@ A kT

∑

where j is summed over all broken bonds at kink site k on edge i. For centrosymmetric growth units, there is just one type of kink site and each kink site is a half-crystal position; so the kink rate, u, is simply Elatt u ¼ R ν ν0 exp ðS 1Þ ð14Þ 2kT Note that eqs 14 and 9 produce eq 6 for the step velocity of centrosymmetric growth units. Also, at equilibrium, S = 1, hence the rate of attachment of the growth units into the kink sites is equal to the rate of detachment of the growth units from the kink sites. For non-centrosymmetric growth units, the expression for kink rate is not linear with supersaturation. It depends on the number of different types of kink sites on the given edge. The kink rate term now accounts for the nonisotropic behavior of the non-centrosymmetric growth units, and hence we need to find the kink incorporation rate to be able to calculate the step velocity. 2791



ARTICLE

Figure 17. Symmetric spiral structure on the {020} face of paracetamol, (a) top view, (b) side view. Asymmetric spiral structure on the {200} face of paracetamol, (c) top view, (d) side view.

Table 3. Asymmetric Spiral Polygon on the {200} Face of Paracetamol edge

critical length (Å)

[001]

1152

[010]

1550

[001]

832

[010]

1550

relative velocity

Table 4. Symmetric Spiral Polygon on the {020} Face of Paracetamol

angle ()

edge

critical length (Å)

relative velocity

angle ()

1.00

1.57

[100]

1491

5.93

1.12

23.71

1.57

[102]

1921

5.92

1.08

152.71

1.57

[001]

1780

1.00

0.93

69.21

1.57

[100]

1491

5.93

1.12

[102]

1921

5.92

1.08

[001]

1780

1.00

0.93

For non-centrosymmetric growth units, we have to deal with more than one type of kink site in series. Let n be the number of different types of kink sites in series in a single orientation f or ). That is, when a growth unit is added to or removed from a site of type Bf, a new kink site of type Af is formed. Conversely, when a growth unit is added to or removed from a site of type Af, a new kink site of type Bf is formed. Similarly, attachments and detachments from A) and B) lead to interconversion of kink

sites between themselves (Figure 13). Therefore, we must correctly account for the transformations of kink sites. The process of incorporation of solute molecules into the kink sites is illustrated in the schematic diagram (Figure 15a). Pk is the probability that the chain is terminated by the growth unit 2792



ARTICLE

Figure 18. (Left) A generalized bonding structure showing double kinks with four different growth units on the face and two different growth units on a single edge; top view. (Right) 3-D view of the same structure.

Table 5. Paracetamol Results no. of face

bond chains

F

lc (Å)

dhkl (Å)

{110}

[110] [001]

0.2 0.014

1325 1780

7.3

[110]

0.2

1325

[001]

0.076

1152

[100]

0.014

1590

[010]

0.21

1370

no yes

{001}

{011}

{201}

{020}

{200}

H-bond H-bond yes no

R

1

1

1

3.05

1

3.22

1

3.55

yes no 6.3

yes

[100]

0.014

1590

[010]

0.21

1370

[100] [122]

0.014 0.3

1640 1870

[100]

0.16

752

[122]

0.3

1870

[102]

0.0014

2730

[010]

0.21

1370

no yes

no 5.3

yes no yes no

5.7

yes

[102]

0.0014

2730

[010]

0.21

1370

[100] [102]

0.172 0.127

1491 1921

[001]

0.014

1780

[001]

0.076

1152

[010]

0.21

1550

no

[001]

0.167

833

no

[010]

0.21

1550

no

no 4.7

yes yes

2

∼20

0

∼20

no 5.8

no

k (that is, the kink site is of the type k and it corresponds to the state k). We want to find an expression for the kink rate on edge i ui ¼ ðRPk νk þ 1 Pk þ 1 Þi

ð15Þ

(corresponding to attachment and detachment) of the solute molecule that is incorporating within the kink site at position k. The first term represents the probability of the site being of type k and the solute molecule incorporating there with the rate R, while the second term represents the probability of the kink site being of type kþ1 (the kink type changes from k to kþ1 upon addition of a solute molecule at site k) and the detachment rate from that particular kink site type, see Figure 15. The kink rate, ui, is the same regardless of which kink sites are chosen in eq 15; that is, the kink rate is for an entire edge, with all its kink sites. Note that, for centrosymmetric growth units, n = 1 so k = 1 as there is just one type of kink site, and the Pk's drop out leading to eq 14. If there are four different types of kinks in series in a single orientation (f or )), the different states go from 1 to 4 in a periodic fashion and then keep repeating. (In other words, k can take values going from 1 to 4; n = 4. The numerical upper bound is still the total number of types of kinks sites, n, in the sense that if the value of the state under consideration k is 4, k þ 1 would be 1. Also if k is 1, then k 1 is 4.) In order to calculate the kink rate using eq 15, we need to calculate the probability Pk of each state k. This is done by setting up a steadystate balance for the influx into and outflux from the state k by addition and removal of solute molecules into and from the different kink sites on edge i. The kth state passes into the (kþ1)st state with rate R and passes into the (k1)st state with rate vk, while the (k1)st state passes into the kth state with rate R and the (kþ1)st state passes into the kth state with rate vkþ1. Figure 15b and eq 16 illustrate this scheme. At steady state, the net rate of transition out of state k is exactly balanced by the net rate into state k. This is represented by a 1D master equation (omitting the index i for clarity) Pk 1 R þ Pk þ 1 νk þ 1 ¼ Pk ðR þ νk Þ

ð16Þ

where (k = 1, 2, 3...n). Shifting indices by 1 in order to get the equation into a more convenient form, we get Pk R þ Pk þ 2 νk þ 2 ¼ Pk þ 1 ðR þ νk þ 1 Þ

The above equation represents the two opposing fluxes 2793

ð17Þ



ARTICLE

Figure 19. The predicted steady-state morphology of paracetamol using our model for non-centrosymmetric growth units at low supersaturation (4%) grown from (a) water, (b) vapor.

Figure 20. The experimental steady-state morphology of paracetamol as observed by Ristic et al.30 (a) Grown from water at low supersaturation; (b) grown from vapor at low supersaturation.

Solving for all the Pk’s, eq 15 gives us a general expression for kink rate on edge i in terms of the attachment rate and detachment rates from different kink sites. ðnÞ

ðRn νi Þ

ui ¼

n

∑

r¼1

ðr 1Þ R n r νi

ð20Þ

where ð0Þ

νi

ð1Þ νi

Figure 21. Experimental growth rates of different faces of paracetamol at different supersaturations measured by Ristic et al.30

Rearranging the equation, R þ νk þ 1 R Pk þ 2 ¼ Pk þ 1 Pk νk þ 2 νk þ 2

ð18Þ

We can generate a system of four equations by inserting different values of k ranging from 1 to 4. Three of the four equations will be independent and we use the condition that n

∑ Pk ¼ 1

k¼1

as our fourth equation.

ð19Þ

ð2Þ νi

ð3Þ νi

¼n

∑ νk k¼1

¼ ¼ ¼

ðr 1Þ νi

!

n

¼ ðν1 þ ν2 þ ν3 þ :::νn Þi i

n

∑ νk νk þ 1

k¼1

! ¼ ðν1 ν2 þ ν2 ν3 þ :::νn ν1 Þi i

n

!

∑ νk νk þ 1 νk þ 2

k¼1

¼

¼ ðν1 ν2 ν3 þ ν2 ν3 ν4 þ :::νn ν1 ν2 Þi i

n

∑ νkνk þ 1 :::νk þ r 2

k¼1

! ð21Þ i

By inserting values for n, we can generate the kink rate expression for any number of different kink types in series. The maximum number of different kink sites theoretically possible on a given edge in a single direction is equal to twice the number of different asymmetric growth units in the unit cell. (The edges are periodic bond chains and hence cell symmetry is repeated after a certain 2794



ARTICLE

Figure 22. (a) Relative growth rates of different faces of paracetamol calculated using the attachment energy model and the steady-state morphology using these relative growth rates (view perpendicular to the c-axis). (b) Steady-state morphology predicted by our new model (view shown perpendicular to the c-axis for comparison).

period; the maximum period can be the number of asymmetric growth units in the unit cell. As far as kink sites are concerned, each growth unit with a given symmetry operator can attach in either of the two orientations, f or ), hence leading to twice the number of kink site types.) n¼2 ui ¼

ðR2 ν1 ν2 Þ ð2R þ ν1 þ ν2 Þ

n¼3 ui ¼

ðR3 ν1 ν2 ν3 Þ ½3R2 þ Rðν1 þ ν2 þ ν3 Þ þ ν1 ν2 þ ν2 ν3 þ ν3 ν1

n¼4 ui ¼ ½ðR4 ν1 ν2 ν3 ν4 Þ=½ð4R3 þ R2 ðν1 þ ν2 þ ν3 þ ν4 Þ þ Rðν1 ν2 þ ν2 ν3 þ ν3 ν4 þ ν4 ν1 Þ þ ðν1 ν2 ν3 þ ν2 ν3 ν4 þ ν3 ν4 ν1 þ ν4 ν1 ν2 Þ

ð22Þ

Knowing the intermolecular bond energies and bond directions, the detachment rates for all the different kink configurations can be calculated. 5.3. Discussion of Unstable Edges. Figure 16 depicts the top view of the {200} face of paracetamol. The specific edge under consideration is the [001] edge. The step front corresponding to this edge is moving in the direction of the arrow. The intermolecular bonding structure on this face is such that the [001] chain is alternating in its edge energy exposed. Thus, this chain has two configurations. To understand the consequences for crystal growth, we first introduce some terminology. The energy perpendicular to the chain, in the direction of motion and falling within the slice is called the edge energy; The energy perpendicular to the chain, in the direction opposite to the direction of motion and falling within the slice is called the reverse edge energy; The energy perpendicular to the chain, and also pointing perpendicular to the face into the solution phase (away from the body of crystal) is called the terrace energy.

The energy perpendicular to the chain and also pointing perpendicular to the face into the body of crystal is called the reverse terrace energy; All energies are local with respect to the edge under consideration. The total energy that then points toward the solution away from the body of the crystal is the sum of the edge energy and the terrace energy, whereas the total energy pointing purely toward the solid state away from the solution phase is the sum of the reverse edge energy and the reverse terrace energy. A stable edge is one that is bound to the crystal with more energy than to the environment, and it will have less tendency to disintegrate compared to an unstable edge in which more energy is exposed to the solution than the body of the crystal. For a stable edge, the net energy required to form kink sites by rearrangement is positive and the net kink rate is positive. For an unstable edge, either the net energy to form kinks by rearrangement from the straight edge is negative implying that the edge would instantaneously rearrange to be completely kinked, or the unstable edge exhibits negative kink rate (dissolution behavior) under supersaturated conditions. Hence, an unstable edge is only a virtual edge. The stable edges would then grow in such a way that only stable edges are formed again, and thus the growth of stable edges occurs in steps of two blocks/monolayers to retain the stable configuration. This is done by forming kinks of double length, as shown in Figure 18. For the particular example of the {200} face of paracetamol, for both configurations of the [001] chain, the energy that is pointing downward from this plane (reverse terrace energy) equals the energy pointing upward from the plane into the solution (terrace energy). Thus, the two [001] edge configurations differ only in their edge energies and reverse edge energies. We label the two configurations as [001]S and [001]U, such that for the [001]S chain, the reverse edge energy is greater than the edge energy making this chain a stable edge. An interesting general observation about a face containing both stable and unstable edges generated due to non-centrosymmetric growth units is that in the polygonized spiral train the opposing edges (i.e., symmetrically related edges whose outward normals point in opposite directions) may not have identical 2795



ARTICLE

Figure 23. (a) Predicted growth rates of different faces of paracetamol from aqueous solvent using our model as a function of the supersaturation ratio in the low supersaturation regime. (b) Relative growth rates of different faces of paracetamol subject to the same conditions. Different colored dots represent pointwise evaluations of the model, and lines represent the best fit polynomial.

Table 6. Lovastatin Solid-State Interactions bond label

bond distance (Å)

dispersive energy (kcal/mol)

coulombic energy (kcal/mol)

total bond energy (kcal/mol)

1

5.97

6.63

0.39

7.02

2 3

11.15 11.39

1.31 1.06

2.15 2.17

3.45 3.23 2.66

4

10.09

2.35

0.31

5

10.38

2.69

0.22

2.47

6

11.19

0.91

0.29

1.19

7

12.01

0.87

0.04

0.91

speeds, thus forming an asymmetric spiral polygon like the one shown in Figure 17c,d. (The relative speeds, critical lengths and angles of the edges are given in Table 3.) This is another noteworthy point of distinction between centrosymmetric and non-centrosymmetric growth units. The faces of a crystal of a centrosymmetric growth unit will only have stable edges in all periodic bond chain directions and will always form a symmetric spiral polygon with opposing edges having equal speeds, like the spiral shown in Figure 17a,b (the{020} face of paracetamol has symmetric bonding structure in the plane and forms a symmetric spiral; see Table 4 for relative speeds, angles, and critical lengths). 5.3.1. Treatment of Unstable Edges. Consider the generalized bonding structure in Figure 18 in which the edge that is built of growth units of type 1 and 3 is the unstable edge configuration

(labeled as u) and the edge consisting of growth units of type 2 and 4 is the stable edge (labeled as s). Since the chains that form alternate stable and unstable edges advance by the formation of double kinks, the master equation (eq 18) developed for the kink rate in Section 5.2.2 is not applicable directly and needs to be modified for these kinds of special edges. (A single edge in Figure 18 has only two growth units in series, whereas if we consider double kinks, we need to consider four different growth units that fill up the kink sites.) The different growth units at the kink sites and the series of events by which they attach to fill up the kink sites is also shown in Figure 18. The growth units attaching to the positions 1 and 3 have very high detachment rates (part of an unstable edge) compared to growth units attaching to positions 2 and 4. So, after the event I in which a growth unit attaches to the 2796



ARTICLE

kink site in position 1, there are two kink sites formed (where the solute growth unit can attach, positions 2 and 3). Since the growth unit has a very low probability of remaining in the kink site at position 3 compared to the position 2, due to the high detachment rate, we assume that the series of events is II, followed by III and then IV. Writing a steady-state balance between the events (note that k takes the values of the event number ranging from I to IV), we obtain a system of equations Pk 1 R þ Pk þ 1 νk þ 1 ¼ Pk ðR þ νk Þ ðk ¼ I; II; III; IVÞ

ð23Þ

which is of the same form as eq 16 with k = 14 (which represented the kink rate involving four different kink sites in series on a single edge) and results in an expression for the kink rate that is same as eq 22 with n = 4.

5.4. Critical Length of Spiral Edges. Predicting the growth rate of a particular face of a crystal requires the determination of the critical lengths for all the spiral edges along with the spiral edge velocities. In the above sections, we discussed the procedure to calculate step velocity. To calculate the critical lengths, we used the approach of carrying out a free energy balance on a single spiral edge, originally proposed by DeYoreo and co-workers47 and further extended by Snyder and Doherty9 and Lovette and Doherty.48 This approach uses the thermodynamic condition that adding a row of molecules of length li greater than the critical length lc decreases the free energy of the system, whereas a row less than the critical length increases the free energy of the system, thereby 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 equation48

lc, i ¼

2ae, i Æφkink, i æ Δμ

ð24Þ

lc, i ¼

2ae, i Æφkink, i æ RT lnðSÞ

ð25Þ

which then becomes

Figure 24. Unit cell structure of lovastatin.

where Δμ is the difference between the chemical potential of the molecules in solution and in the crystal, S is the supersaturation ratio, ae,i is the distance between molecules along edge i, and is the arithmetic average kink energy calculated

Figure 25. Bonding structure of the relevant F-faces of lovastatin (a) {210} face, (b) {020} face, (c) {011} face, (d) {101} face. 2797



ARTICLE

Figure 26. Predicted steady-state morphology of lovastatin.

Figure 27. The experimental steady-state morphology of lovastatin grown from (a) isopropanol, (b) methanol.

velocities. All these quantities thus depend on the nature of the solvent and this directly impacts the final steady-state shape of the crystal. We scale the bonds according to the approach described by Winn and Doherty29,49 and extended by Snyder and Doherty.9 The method is a classical one based on estimating solid and solvent interactions independently (making use of solubility parameters for the solvent phase50) and calculating the interfacial energy from the adhesive and cohesive components.

6. RESULTS

Figure 28. The attachment energy morphology of lovastatin.

according to the procedure in Section 5.2.1. Critical lengths are expected to be in the range from 10 to 1000 nm. 5.5. Solvent Effect. While the crystal lattice is growing, solute molecules that are surrounded by the fluid phase (solvent) are preferentially incorporating within the kink sites. Since the solute phase is constantly surrounded by the solvent, we cannot use the intermolecular bond strengths calculated by the AMBER force field directly. All the solidsolid broken bonds need to be interfaced to account for the effect of the solvent. These broken bond strength values are used in the calculation of edge energies, kink energies, kink rates, and hence the critical lengths and step

6.1. Steady-State Morphology of Paracetamol. Using the methods and equations described in the earlier sections, the step velocities and critical lengths of the spiral edges on each relevant face of monoclinic paracetamol (P21/a space group) that can grow layer-by-layer via spirals were estimated, and hence it was possible to calculate the relative growth rates of faces. Some of the absolute quantities that can be predicted are presented in Table 5. (For the {020} faces, the opposing edges, that is, the outward normals pointing in opposite directions and related by symmetry, have exactly the same values and so only the positive vector directions are reported.) The detailed results for all of the intermediate data (kink energies, edge energies, kink rates, step velocities, etc.) can be found in the Supporting Information. Upon the basis of the predicted relative growth rates of the faces, the steady-state morphologies of paracetamol grown from water and from the vapor at a supersaturation of 4% were constructed and are shown in Figure 19a,b. The predicted shapes are compared to the corresponding experimental shapes measured 2798



ARTICLE

Figure 29. Bonding structure for the {210} face of lovastatin.

Table 7. Configuration of Stable and Unstable Edges of the [001] Chain on the {210} Face of Lovastatina type and order of edge

ER

TR

E

T

Table 8. Configuration of Stable and Unstable Edges of the [001] Chain on the {210} Face of Lovastatin

ER þ TR E þ T S or U

1 2

5.92 8.02 6.46 7.42 6.46 7.42 5.92 8.02

13.94 13.88

4

5.92 8.02 5.32 8.56

3

5.32 8.56 5.92 8.02

type and order of edge ER 1 3

TR

E

T

13.88 13.94

S U

6.46 8.02 5.92 7.42 5.92 8.56 5.322 8.02

13.94

13.88

S

4

5.32 8.02 5.92

8.56

13.88

13.94

U

2

5.92 7.42 6.46

8.02

ER þ TR E þ T S or U 13.34 13.34

S S

13.34

14.48

U

13.34

14.48

U

14.48 14.48

a

ER stands for reverse edge energy, TR is the reverse terrace energy, E is the edge energy, and T the terrace energy (all absolute values). S and U denote whether the given chain configuration is a stable or an unstable edge, respectively.

by Ristic and co-workers30(Figure 20). We see that the predicted morphology of paracetamol at low supersaturation (4%), grown from water and vapor, matches remarkably well to the experimental shapes. Also, as one can see from Table 5, the {110} face is the slowest growing face, followed by {001}, and {201} faces, respectively, which is the same trend observed by Ristic and coworkers for aqueous growth at low supersaturations; see Figure 21. Also, note that the predictions of the attachment energy model (shown in Figure 22a) gives a different morphology with multiple facets (compare with Figure 22b), thus demonstrating that a more rigorous and mechanistic model than the attachment energy model is needed to correctly capture the crystal growth behavior. We see from Figure 21 that the relative growth rate of the {110} face is very low compared to the other faces in the low supersaturation range (of 010%). A school of thought exists that this face exhibits a dead zone at low supersaturation with recovery at higher supersaturation due to selective absorption of impurities on this face.51 However, our analysis is that the {110} face is neither uniquely hydrophobic nor uniquely hydrophilic to behave very differently from other faces. Our explanation is that the relative growth rate of the {110} face is significantly low compared to the other faces at low supersaturation primarily due to its bonding structure (one of the periodic bond chains on this face, the [110] chain, has eight different types of kink sites). This causes the [110] chain to have very low kink rates and hence a very low step velocity and consequently a very low relative growth rate of the {110} face. The growth rate of this {110} face becomes higher at higher supersaturation (with a corresponding decrease in the morphological importance).30 This sudden increase in growth rate with supersaturation and a different morphology is believed to occur after a critical value of supersaturation due to a change in

mechanism from spiral growth to 2-D nucleation. The proof of concept of this hypothesis of mechanism change is the subject of a future paper.27 The model can also predict the growth rates (to within a constant multiple) at different supersaturations, as shown in Figure 23a,b. Comparing the predicted growth rates with those measured by Ristic (Figure 21), we see that the scales agree in the low supersaturation regime (0 to 10%) to within a constant multiple; the growth rates increase from 0 to 1. The morphology remains the same in the low supersaturation regime, as reflected by the constant relative growth rates of different faces (Figure 23b), in line with the experimental observations of Ristic,30 and of Shekunov and Grant.31 The actual experimental values for G are uncertain because different groups report different growth rates,30,31 reflecting the fact that they are difficult to measure. Hence, based on the above analysis, the model is compatible with what is known from experiments. Regression analysis of the predicted growth rate of faces vs supersaturation indicates that growth rate of the face, G0 , scales almost exactly as (S 1)2. Thus, the relative growth rates are expected to be independent of supersaturation since, Ri = (Gi0 /G0 110) = (ai(S 1)2)/(a110(S 1)2) = (ai/a110). It appears that the relative growth rates remain constant over the range of supersaturations where the spiral growth model applies. We expect that changes in relative growth rates are an indication of a change in growth mechanism. 6.2. Morphology of Lovastatin. We also tested our spiral growth model on lovastatin, a higher molecular weight pharmaceutical compound (molecular weight = 404.5) with a more complicated chemical bonding structure than paracetamol. Lovastatin is a member of the drug class of statins, used for lowering cholesterol (hypolipidemic agent) in those with hypercholesterolemia and so for preventing cardiovascular disease. Lovastatin crystallizes in the space group P212121 with the lattice parameters a = 22.15 Å, b = 17.32 Å, c = 5.97 Å, and R = β = γ = 90 . There are four molecules per unit cell having different symmetry operators, shown in Figure 24 by different colors. 2799



ARTICLE

Figure 30. Stable and Unstable edge configurations for the advance of the [001] edge on the {210} face of lovastatin.

The calculations were done using the same procedure as for paracetamol and the PBCs and the F-faces were calculated. The relevant F-faces for this system are {210}, {101}, {011}, and {020} faces. The calculated lattice energy for lovastatin using the AMBER force field is 44.0 kcal/mol. Lovastatin exhibits a different configuration of stable and unstable edges that is an extension of the case discussed in Section 5.3.1. On the {210} face, the crystal structure repeats after four monolayers of the [001] step. The bonding structure on this face is such that there are two stable [001] edge configurations together followed by two unstable [001] edge configurations, and this goes on repeating. This case is dealt with in Appendix A, and the details about bond chains parallel to each of the above F-faces and relative growth rates are provided in the Supporting Information. Using these relative growth rates, the steady-state morphology of lovastatin is predicted and is shown in Figure 26. Lovastatin was predicted to grow as long needles from alcohols. Since no detailed experimental morphology is available in the literature, we performed our own experiments of growing lovastatin in two different solvents, isopropanol and methanol, by the cooling crystallization method. The crystals were grown at low supersaturation (510%) at 16 C in a small, quiescent crystallizer equipped with a video microscope, moveable stage, and Peltier cell for precise temperature control (experimental setup is described elsewhere52,53). We wish to remark here that the experiments were started only after the predictions were completed. Thus, the prediction of the morphology of lovastatin was made without any prior knowledge of the experimental shape, thus providing a rigorous test of the predictive power of the model. As seen in Figure 27, lovastatin grows as very long needles in both the solvents of different polarities, which is in line with the predicted morphology. We again compare the experimental morphology with the attachment energy model, as shown in Figure 28, and observe that the attachment energy model does not predict needles.

7. CONCLUSIONS We have developed a mechanistic model that can predict the crystal morphology of any molecular organic crystal, centrosymmetric or non-centrosymmetric, given the crystal structure of the polymorph of interest. As this model is based on the actual mechanism of spiral growth, it is a more natural and accurate means of predicting crystal shape. It is also much faster to execute than molecular simulations. At the same time, there are challenges for carrying out such mechanistic modeling predictions. Input system preparation is an important starting point for morphological predictions. Many researchers have studied the sensitivity of the selected forcefield to the accuracy of morphology prediction for organic compounds,54 including paracetamol.33 The highest fidelity test for selecting a particular force field is to compare rms deviations of the structural lattice parameters as well as the energy differences, by lattice energy minimization methods, as described by Brunsteiner et al.54 and Wang et al.55 A more direct and easier method of choosing the force field is to compare the calculated lattice energy to the experimental value. We use this direct approach wherever possible. Differences between the experimental lattice energy and the heat of sublimation of less than about 10 kJ/mol are not significant,56 as these are likely to represent thermodynamic approximations and experimental error. We have used the AMBER force field for both paracetamol and lovastatin, as this force field was specifically parametrized for organic systems and is known to work very well with very low rms deviations in the structure for drug molecules.55 We have argued that the BCF assumption of a Kossel crystal is to be used with caution as it is applicable only for simple centrosymmetric molecules, whereas for many other systems of realistic complexity, we must account for the Non-Kossel behavior. We have developed generalized expressions for kink density and kink incorporation rates for any number of different types of 2800


Crystal Growth & Design growth units and kink sites, which can be utilized for the calculation of step velocity. The predictive capability of the model has been demonstrated using examples of two industrially relevant systems, paracetamol and lovastatin, where we obtained very good agreement between the predicted and experimental steady-state morphologies. We remark that the general model is not limited to just one chemical entity within the unit cell (molecular crystals). Co-crystals have two different kinds of molecules (chemically different) existing in specific sites within the same crystal lattice. Designing co-crystals is the research of the future for pharmaceutical companies as it allows them to combine multiple drugs into a single crystal, making it more effective. Common cocrystals include hydrates and solvates where water or another solvent are integrated into specific positions in the lattice. Cocrystals are necessarily non-Kossel crystals with multiple kink sites, with the additional constraint of having chemically different growth units within a crystal. For predicting crystal morphology of salts, solvates and cocrystals, we need to make appropriate modifications to the force field such that the intermolecular bonding is correctly captured. This is a subject of ongoing research.

’ APPENDIX A Unstable Edge Configuration on the {210} Face of Lovastain. Figure 29 shows the bonding structure for the

{210} face of lovastatin. The direction of the arrow signifies the direction of advance of the [001] edge, whereas the [001] edge advances in the direction opposite to the arrow. As seen in the diagram, there are four chain configurations for each of the [001] and [001] edges (1, 2, 3 and 4, shown by the red colored bond chains and each has a same bond strength of 7.02 kcal/ mol). The stable and unstable edge configurations are computed in Tables 7 and 8 by calculating which edge configurations have greater energy pointing within the body of the crystal than to the environment. As seen from Table 7, the [001] chain has alternate stable and unstable configurations, a case identical to which has already been described in the main text (Section 5.3.1). However, the [001] chain configurations (Table 8) advance such that there are two stable edges followed by two unstable edges. We now find an expression for the kink rate for this case. Consider Figure 30a. Even though each of the edges are comprised of different growth units, the total energy pointing towards the body of crystal (ER þ TR) and the total energy directed towards the solvent (E þ T) is the same for edges 2 and 4 (also 1 equals 3). Hence, the four edge configurations simplify to just two different edge configurations, unstable and stable. (The unstable configuration is shown by the red blocks in Figure 30b,c and has a higher rate of detachment, meaning it has less tendency to join the crystal body compared to the stable edge, shown by the green blocks in those figures.) Thus, Figure 30a is equivalent to the simplified form Figure 30b. The individual sequence of steps by which the kink sites would be filled up is shown in Figure 30c. Event I shows filling up of the single kink site by the red growth unit which generates two kink sites that are filled up by two red growth units as shown in Event II. This generates three kink sites; however, the detachment rate of the red growth unit is much higher than the green growth unit. This means that the green growth unit has a significantly greater tendency of attaching than the red growth units, so we assume a sequential attachment of green growth unit

ARTICLE

followed by the red growth units as shown in the subsequent events. Also, Events IV, V, and VI are identical to I, II, and III, respectively, so there are just three different states. Writing a steady-state balance between these states, similar to the approach used in Section 5.3.1, we get PII R þ PII ðν2 þ ν2 Þ ¼ PI ðR þ RÞ þ PIII ν1 PIII R þ PIII ν1 ¼ PII R þ PI ν1 PI ðR þ RÞ þ PI ν1 ¼ PIII R þ PII ðν2 þ ν2 Þ

ð26Þ

where v1 and v2 signify the detachment rates of the growth units comprising the unstable edge (red) and the stable edge (green), respectively. We obtain an expression for the kink rate on the [001] edge by defining it as follows u ¼ PIII R PI ν1

ð27Þ

Solving eq 26 and making use of the Pk's in eq 27, we obtain u¼

ð2R3 2ν21 ν2 Þ ½5R2 þ Rð3ν1 þ 2ν2 Þ þ ν21 þ 4ν1 ν2

ð28Þ

’ ASSOCIATED CONTENT

bS

Supporting Information. Tables 1 and 2 list intermediate values used for constructing steady-state morphologies of paracetamol and lovastatin. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*Phone: 805-893-5309; fax: 805-893-4731; e-mail: mfd@ engineering.ucsb.edu.

’ ACKNOWLEDGMENT We are grateful for the financial support provided by Abbott Laboratories and Eli Lilly. Dr. N. Variankaval of Merck & Co. graciously provided samples of lovastatin as well as its crystal structure. ’ REFERENCES (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) Rimer, J. D.; An, Z.; Zhu, Z.; Lee, M. H.; Goldfarb, D. S.; Wesson, J. A.; Ward, M. D. Science 2010, 330, 337–341. (4) Yip, C. M.; Ward, M. D. Biophys. J. 1996, 71, 1071–1078. (5) Sours, R.; Zellelow, A.; Swift, J. J. Phys. Chem. B 2005, 109, 9989–9995. (6) Vekilov, P.; Alexander, J. Chem. Rev. 2000, 100, 2061–2090. (7) Burton, W. K.; Cabrera, N.; Frank, F. C. Phil. Trans. R. Soc. A 1951, 243, 299–358. (8) Chernov, A. A. Modern Crystallography III. Crystal Growth; Springer-Verlag: Berlin, 1984. (9) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. A 2009, 465, 1145–1171. (10) Zhang, J.; Nancollas, G. H. J. Colloid Interface Sci. 1998, 200, 131–145. 2801


Crystal Growth & Design (11) Chernov, A.; Rashkovich, L.; Vekilov, P. J. Cryst. Growth 2005, 275, 1–18. (12) Rashkovich, L.; Petrova, E.; Chernevich, T.; Shustin, O.; Chernov, A. Crystallogr. Rep. 2005, 50, S78–S81. (13) Cuppen, H.; Meekes, H.; van Enckevort, W.; Vlieg, E. Surf. Sci. 2004, 571, 41–62. (14) Cuppen, H.; Meekes, H.; van Enckevort, W.; Vlieg, E. J. Cryst. Growth 2006, 286, 188–196. (15) Kaischew, R.; Budevski, E. Contemp. Phys. 1967, 8, 489–516. (16) Voronkov, V. V. Sov. Phys. Cryst. 1973, 18, 19–23. (17) Teng, H. H.; Dove, P. M.; Orme, C. A.; De Yoreo, J. J. Science 1998, 282, 724–727. (18) Case, D. A. AMBER 10; University of California: San Francisco, 2008. (19) 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. (20) Hartman, P.; Perdok, W. G. Acta Crystallogr. 1955, 8, 49–52. (21) Frank, F. C. In Growth and Perfection of Crystals; Doremus, R. H., Roberts, B. W., Turnbull, D., Eds.; Wiley: New York, 1958. (22) Chernov, A. A. Sov. Phys. Cryst. 1963, 7, 728–730. (23) Zhang, Y.; Sizemore, J. P.; Doherty, M. F. AIChE J. 2006, 52, 1906–1915. (24) Kossel, W. Nachur. Ges. Gottingen 1927, 206, 135–145. (25) Markov, I. V. Crystal Growth for Beginners, Fundamentals of Nucleation, Crystal Growth and Epitaxy; World Scientific: Singapore, 2003. (26) Herbstein, F. H.; Schoening, F. R. L. Acta Crystallogr. 1957, 10, 657–663. (27) Lovette, M. A.; Doherty, M. F. 2011 (in preparation). (28) Sizemore, J. P.; Doherty, M. F. J. Cryst. Growth 2010, 312, 785–792. (29) Winn, D.; Doherty, M. F. AIChE J. 1998, 44, 2501–2514. (30) Ristic, R.; Finnie, S.; Sheen, D.; Sherwood, J. J. Phys. Chem. B 2001, 105, 9057–9066. (31) Shekunov, B.; Grant, D. J. Phys. Chem. B 1997, 101, 3973–3979. (32) Boerrigter, S. X. M.; Cuppen, H. M.; Ristic, R. I.; Sherwood, J. N.; Bennema, P.; Meekes, H. Cryst. Growth Des. 2002, 2, 357–361. (33) Cuppen, H.; Day, G.; Verwer, P.; Meekes, H. Cryst. Growth Des. 2004, 4, 1341–1349. (34) Lu, J. J.; Ulrich, J. Cryst. Res. Technol. 2005, 40, 839–846. (35) Beyer, T.; Day, G. M.; Price, S. L. J. Am. Chem. Soc. 2001, 123, 5086–5094. (36) Haisa, M.; Kashino, S.; Kawai, R.; Maeda, H. Acta Crystallogr. 1976, B32, 1283–1285. (37) Clydesdale, G.; Docherty, R.; Roberts, K. J. Comput. Phys. Commun. 1991, 64, 311–328. (38) Lovette, M. A. Doctoral Dissertation, University of California Santa Barbara, 2011 (in preparation). (39) Frisch, M. J. Gaussian 03, Revision C.02. 2003 (40) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269–10280. (41) Li, T.; Feng, S. Pharm. Res. 2006, 23, 2326–2332. (42) Hartman, P.; Perdok, W. G. Acta Crystallogr. 1955, 8, 521–524. (43) Grimbergen, R.; Reedijk, M.; Meekes, H.; Bennema, P. J. Phys. Chem. B 1998, 102, 2646–2653. (44) Chernov, A. A. J. Cryst. Growth 2004, 264, 499–518. (45) Choi, C. S.; Prince, E. Acta Crystallogr. Sect. B 1972, 28, 2857–2862. (46) Gilmer, G. H.; Bennema, P. J. Appl. Phys. 1972, 43, 1347–1360. (47) Thomas, T. N.; Land, T. A.; Martin, T.; Casey, W. H.; DeYoreo, J. J. J. Cryst. Growth 2004, 260, 566–579. (48) Lovette, M. A.; Doherty, M. F. J. Cryst. Growth 2011, accepted for publication. (49) Winn, D.; Doherty, M. F. AIChE J. 2000, 46, 1348–1367. (50) Barton, A. F. M. Chem. Rev. 1975, 75, 731–753. (51) Ristic, R. I.; DeYoreo, J. J.; Chew, C. M. Cryst. Growth Des. 2008, 8, 1119–1122.

ARTICLE

(52) Veesler, S.; Lafferrere, L.; Garcia, E.; Hoff, C. Org. Process Res. Dev. 2003, 7, 983–989. (53) Lovette, M. A.; Muratore, M.; Doherty, M. F. AIChE J. 2011, accepted for publication. (54) Brunsteiner, M.; Price, S. L. Cryst. Growth Des. 2001, 1, 447–453. (55) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157–1174. (56) Gavezzotti, A.; Filippini, G. In Theoretical Aspects and Computer Modeling. the Molecular Solid State; Gavezzotti, A., Ed.; John Wiley & Sons Ltd.: New York, 1997; Chapter Energetic aspects of crystal packing: Experiment and computer simulations, pp 6297.

2802


Spiral Growth Model for Faceted Crystals of Non ... - ACS Publications

Recommend Documents