## Abstract

We present a systematic modelling methodology using the spiral growth mechanism of Burton, Cabrera and Frank to predict the steady-state shape of organic molecular crystals grown from solution. This methodology has been developed to eliminate the need for special modifications for each new crystal system studied. Therefore, the mechanisms and choices for spiral shapes, edges and evolution are mathematically determined as governed by the underlying solid-state chemistry and physics. The power of the approach is demonstrated for several crystal systems: naphthalene grown from both ethanol and cyclohexane; anthracene grown from 2-propanol; and glycine grown from water. The predicted crystal shapes are in good agreement with experiment.

## 1. Introduction

It has been nearly 60 years since Burton *et al*. (1951) first published their landmark work, predicting the growth of crystals via a spiral mechanism emanating from a screw dislocation. Since that time, many detailed experiments have observed this mechanism in action for a wide range of crystallizing materials (Yip & Ward 1996; Paloczi *et al*. 1998; Vekilov & Alexander 2000; Sours *et al*. 2005; Ranguelov *et al*. 2006). As a result, the Burton, Cabrera and Frank (BCF) spiral growth mechanism is recognized as the most important modelling tool in the hands of crystal growers. The model itself has undergone improvements over the years; nonetheless, its implementation remains essentially an art owing to the large number of intuitive decisions that must be made, in particular for faceted crystals. In this paper, we report a systematic procedure for implementing the BCF spiral growth mechanism to predict crystal morphology.

Prediction of crystal shape continues to be an important topic of scientific interest owing to the impact that crystal shapes have on the properties of materials, including dyes, catalysts, semiconductors and pharmaceuticals. Crystal shape affects both the performance and processing of the material. Different crystal faces can have different chemical properties such as surface reactivity or hydrophobicity, and, since the crystal shape determines the relative area of each face, it also determines the magnitude of these properties. Additionally, bulk dissolution rate and mechanical properties are dependent on crystal shape. The crystal shape can also impact downstream processes such as filtering, washing or drying, as well as slurry viscosity and bulk density. Because the desired shape is often specific to the industry and the application for each particular crystal system, a general methodology for predicting shape as a function of the important design parameters is desired.

The first models that were used to predict crystal growth shapes can be traced to the Bravais–Friedel–Donnay–Harker model (Bravais 1886; Fridel 1907; Donnay & Harker 1937), where the only input for crystal shape prediction is the interplanar spacing. Later, Hartman & Perdok (1955) developed the attachment energy model to explicitly include the energetic interactions required to attach a new slice of molecules to a surface. In the past decade, the state-of-the-art methods in modelling have progressed to include the actual mechanisms for growth. Simulation techniques include the Monte Carlo algorithm MONTY (Boerrigter *et al*. 2004) and the kinetic Monte Carlo studies of Piana & Gale (2005) and Gilmer (Zepeda-Ruiza *et al*. 2006). Additionally, a class of mechanistic kinetic models, based upon the methods of Burton *et al*. (1951) and Chernov (1984), was developed by Winn & Doherty (1998). These kinetic models incorporate the effect of solvent, as well as the anisotropic energetic interactions within crystals (Heng *et al*. 2006), and have the potential to be implemented rapidly, making them uniquely suited to both product and process development. Despite these advances in first-principles modelling, there are still significant challenges to overcome before they will be widely used.

In this paper, we report the next-generation spiral growth model for faceted crystals. The model significantly reduces the need for special modifications to be made for each crystal system, thus unifying the theory. Unification is achieved by strengthening the connection between spiral motion and the underlying solid-state chemistry. The model has been automated, which enables future consideration of more complicated systems such as complex pharmaceutical molecules, co-crystals and the effect of mixed solvents and additives.

New methodologies are reported in several areas. The model specifically accounts for crystal geometry (distance of growth propagation, detailed edge and kink area calculations and step heights) along with the physical basis for its inclusion. The inclusion of important intermolecular interactions is based solely upon energetics rather than distances, improving on previous formulations (Winn & Doherty 1998). Also, the selection of step edges to include in the spiral growth rate calculation is mathematically determined. The critical lengths of edges on spirals are now determined using a model based on the open geometry of a spiral rather than the classical Gibbs–Thomson analysis of the equivalent closed two-dimensional shape. Not only is each of these new methodologies incorporated into the model and algorithm in order to reduce the number of ad hoc decisions, but is also applied uniformly across each crystal system, thus providing for systematic predictions. Despite this added modelling fidelity, the full calculation can still be performed rapidly (tens of seconds to a few minutes) using an automated Matlab program on a desktop PC with a single processor that takes crystallographic and molecular interaction input data and generates the resulting steady-state crystal shape, as well as all intermediate data.

This paper first describes the modelling theory with an emphasis on new aspects of spiral dynamics. Then, two example systems, which have not previously been studied with these kinetic models, are reported and compared with the experiment: naphthalene grown from cyclohexane and ethanol, and anthracene grown from 2-propanol. Finally, we perform a calculation for α-glycine grown from water. This calculation was previously performed (Bisker-Leib & Doherty 2003) with many of the solid-state chemistry decisions argued intuitively rather than systematically determined.

## 2. Growth spirals

For crystals growing by the steady flow of steps across surfaces, the mathematical expression for the rate of growth normal to a crystal surface (*hkl*) is(2.1)where *v*_{i} is the perpendicular step velocity across the surface in the *i*th edge direction; *h* is the step height; and *y*_{i} is the interstep distance. Each of these quantities is face dependent; however, the subscript *hkl* will be dropped for the ease of presentation. Prediction of the growth rate of a face growing by the spiral mechanism thus requires the independent determination of each of these quantities for each crystal face. One alternative expression for equation (2.1), which is often convenient for interpreting experiments (Sours *et al*. 2005), is(2.2)where *ρ*_{i} is the slope of the spiral hill averaged over many steps. It should be noted that the normal growth rate of a given face is independent of the edge direction selected in equations (2.1) or (2.2). (While *v*_{i} and *y*_{i} are dependent on the edge, their ratio is not; thus, the growth rate is the same for each value of *i*.) Equation (2.1) can be reformulated as(2.3)where *τ* is the characteristic rotation time of the spiral. This characteristic time corresponds to the time between consecutive passes of a step across a point on the face. The characteristic spiral rotation time can be calculated by the time that it takes for each edge to appear in one full first turn of the spiral. Beyond the first turn, the time it takes for yet another pass of a step to cross the same point on the face is determined solely by the creation of another first turn. Thus, the rotation time is only relevant during the first turn. In some complex systems, consecutive layers of a crystal can have alternating properties, including rotation time. In these cases, it is necessary to evaluate a number of turns corresponding to one full grouping of the growth pattern. From a modelling perspective, equation (2.3) is the most important and convenient formulation.

We begin by assuming that the velocity of a step is governed by(2.4)which was first proposed by Kaischew and later by Voronkov (Kaischew & Budevski 1967; Voronkov 1973). This profile indicates that, when a spiral edge is shorter than its critical length, it will not grow at all, and once it is above its critical length, it will grow with a constant step velocity. The concept of a critical length for growth has been experimentally validated for calcite (Teng *et al*. 1998) and lysozyme (Chernov *et al*. 2005). These results also demonstrate that equation (2.4) is a first good approximation for the velocity profile. Voronkov's conditions in equation (2.4) are not reversible; if an edge begins to shrink in length, its perpendicular velocity will continue to be *v*_{∞} as the edge length passes below *l*_{c} and smoothly shrinks to zero.

### (a) Spiral edge determination and first turn dynamics

For a convex spiral, the rotation time is given by(2.5)where *α*_{i,i+1} is the angle between edges *i* and *i*+1; *l*_{c,i+1} is the critical length of edge *i*+1; and *N* is the number of relevant spiral edges. (*N* will be some or all of the total set of *M* edges on a given face. A method for determining which edges in *M* to include in *N* follows.) We define a numbering scheme in which the first edge to grow is labelled 1 and the last is *N*. Therefore, side *i* exposes side *i*+1. This numbering scheme is used for the spirals rotating in either a clockwise or anticlockwise direction. The rotation direction is defined by the advancement of the outward normal of edge *i*+1 relative to that of edge *i*. Equations (2.3) and (2.5) are used to predict the growth rate of a crystal face. Thus, rather than directly predicting the step velocity and interstep distance for a single-step direction, as required for equation (2.1), we predict the critical length and the step velocity for each edge of the spiral. (The interstep distance is then given by *y*_{i}=*v*_{i}*τ*.) Additionally, we must determine the correct number of edges to include in the set *N*. Once *N* is known and the crystallography is given, the angles between these spiral edges, *α*_{i,i+1}, are fixed and easily determined. As shown by equation (2.3) and equation (2.5), the growth rate of a face is determined by the characteristic time for the first turn of the spiral. Thus, we first focus attention on the initial spiral turn, its dynamics and which edges of the spiral are important for the calculation of *τ*. Additionally, the model for the first turn extends to incorporate spirals away from the dislocation and leads to interesting insight into their behaviour.

The shape of the spiral is determined by the relationship between the critical lengths and the velocities of each spiral edge, as well as the angles between them. Given the perpendicular step velocity of each spiral edge and the angles between them, the corresponding tangential velocity (the rate at which the edge lengths increase or decrease) of each spiral edge is also known. The tangential velocity of a spiral edge *i* is given by(2.6)where *α*_{i,j} is the angle between edges *i* and *j*; is the tangential velocity of edge *i*; and *v*_{i} is the perpendicular step velocity of edge *i*. This is in direct analogy with the representation of the edge of a two-dimensional crystal, as was first published by Kozlovskii (1957). While a spiral is not a closed geometric object, in contrast to a two-dimensional crystal, the calculation is only a function of the edge of interest, along with each adjacent edge. Thus, the edge under consideration does not have to be part of a closed object for this equation to be valid.

Since the step velocity profile is assumed to be governed by equation (2.4), each edge of the spiral undergoes three different stages of tangential growth. Assuming that an edge (*i*−1) already exists and its length is greater than its critical length (*l*_{i−1}>*l*_{c,i−1}), these three stages (for edge *i*) correspond to (i) when edge *i* begins to appear and its length is less than its critical length (*v*_{i}=0) and edge *i*+1 has yet to appear (*v*_{i+1}=0) (figure 1*a*), (ii) when edge *i* is larger than its critical length and edge *i*+1 has appeared but is less than its critical length (*v*_{i+1}=0) (figure 1*b*), and (iii) when edge *i* and *i*+1 are both larger than their critical lengths (figure 1*c*).

During the first stage, edge *i* is less than its critical length *l*_{c,i}. Thus, the only non-zero step velocity in equation (2.6) is *v*_{i−1} (edge *i*+1 has yet to appear) and the tangential velocity reduces to(2.7)Once edge *i* is larger than its critical length, it begins to move, exposing edge *i*+1, which is initially shorter than its critical length *l*_{c,i+1}; thus *v*_{i+1}=0. The tangential velocity expression for edge *i* is then given by(2.8)Finally, once edge *i*+1 begins to move, edge *i* obtains its final tangential velocity given by equation (2.6) . During the second stage, the new edge *i*+1 is being exposed and has a step velocity of zero. Thus, if the step velocity of edge *i* is large, then edge *i* could have a negative tangential velocity during this second phase (i.e. it shrinks in length), while it must have had a positive tangential velocity during the first phase (and either a positive or negative tangential velocity in the third stage). If step *i* moves too fast, it is even possible for edge *i* to disappear during stage (ii) before edge *i*+1 has begun to move. If this is the case, the spiral rotation time will not depend on edge *i*, and it should not be included in the calculation of the spiral time in equation (2.5).

In order to mathematically determine whether spiral edge *i* disappears before edge *i*+1 begins to move, a comparison is made between the time it would take for edge *i* to disappear and the time it takes for edge *i*+1 to reach its critical length . In some cases, *t*_{dis}<0, implying edge *i* would never disappear owing to the geometry of the spiral; thus, edge *i* should be included in the set of *N* edges for calculating *τ*. An example where *t*_{dis}<0 is a square spiral where since the geometry (see equation (2.8) with all angles equal to 90°) does not allow for the tangential velocity of any edge to ever be negative. In other cases, edge *i* may eventually disappear (*t*_{dis}>0). For such cases, if the time for edge *i* to disappear is less than the time it takes for edge *i*+1 to reach its critical length (*t*_{dis}<*t*_{crit}), then edge *i* should not be included in the spiral time calculation, otherwise it should be included. By equating these two times and solving for the velocity of edge *i*, a velocity for edge inclusion can be determined. This inclusion velocity for edge *i* of the spiral is given by(2.9)Thus, when the step velocity of edge *i*, *v*_{i}, is greater than or equal to its inclusion velocity, , (and *t*_{dis}>0), it will disappear during phase 2 of its growth and should not be included in the spiral rotation time calculation. Whereas if the step velocity of edge *i* is less than this inclusion velocity (or if *t*_{dis}<0), it should be included in the spiral rotation time calculation. Thus, the number of edges *N* to include in the spiral rotation time equation is equal to the number of edges that pass the inclusion velocity test of equation (2.9), and the crystal growth rate prediction reduces to a prediction of critical lengths and velocities for each spiral edge on each crystal face.

Figure 2 demonstrates the differences between two spirals, each with six edges considered, but one which contains two edges that do not persist beyond the growth of the following edge. Figure 2*a* corresponds to the resulting spiral for the edge velocities (in arbitrary units) *v*_{1}=*v*_{4}=2.0, *v*_{2}=*v*_{5}=1.0 and *v*_{3}=*v*_{6}=1.5, with critical lengths *l*_{c,1}=*l*_{c,4}=2.0, *l*_{c,2}=*l*_{c,5}=3.0 and *l*_{c,3}=*l*_{c,6}=4.0 and angles *α*_{1,2}=*α*_{4,5}=50°, *α*_{2,3}=*α*_{5,6}=60° and *α*_{3,4}=*α*_{6,1}=70°. Figure 2*b* corresponds to the resulting spiral for *v*_{1}=*v*_{4}=2.0, *v*_{2}=*v*_{5}=1.5 and *v*_{3}=*v*_{6}=1.0, with critical lengths *l*_{c,1}=*l*_{c,4}=2.0, *l*_{c,2}=*l*_{c,5}=4.0 and *l*_{c,3}=*l*_{c,6}=3.0 and angles *α*_{1,2}=*α*_{4,5}=70°, *α*_{2,3}=*α*_{5,6}=60° and *α*_{3,4}=*α*_{6,1}=50°. For the spiral in figure 2*a*, each of the edges of the spiral passes the test to determine whether it should be included in the calculation of *τ*. However, for the spiral in figure 2*b*, the velocities of edges 1 and 4 are greater than those of the inclusion velocities (and *t*_{dis}>0); thus, these two edges are not included in the calculation of *τ* for that spiral. The spiral rotation time for the first spiral is then given by equation (2.5) with *N*=6 (*τ*≃11.7), whereas that for the second spiral is given by(2.10)(*τ*≃10.4). The smaller value of the rotation time leads to a larger value of the growth rate, *G*_{i}, for the same face.

### (b) Spiral shape evolution

Based upon equations (2.3) and (2.5), the growth rate of a crystal face is determined purely by the dynamics of the first turn. Once the first turn is completed, the interstep distance in each edge direction between one rotation and the next is fixed. Also, the step velocity is assumed to be constant when an edge is longer than its critical length and the height of a step is constant. According to equation (2.1), the face growth rate does not depend on anything beyond the first turn, since each of the terms in the equation is fixed. Nonetheless, the shape of the spiral evolves as it moves away from the dislocation at its centre, and eventually achieves a steady-state spiral shape. Spiral edges can also disappear during this evolution process after they reach their third stage of growth (where edge *i*+1 has also begun to grow), but it is rather unlikely that new edges would appear without a change in the system (e.g. a change in temperature or the addition of a surface active impurity into the solution). Nevertheless, even the disappearance of spiral edges in (and beyond) their third stage of growth will not affect the normal growth rate of a face.

Since the tangential velocity of a spiral edge is determined by the same expression as the tangential velocity of a two-dimensional crystal, the critical step velocity that determines whether a spiral edge will disappear at some distance away from the centre of the spiral (beyond stage (ii) of its growth) is also given by the same expression. This critical step velocity for the disappearance of a spiral edge beyond its second stage of growth is given by(2.11)which is identical to the expression derived by Szurgot & Prywer (1991) for the two-dimensional crystals. As a spiral evolves in shape and extends outwards from the dislocation, a steady-state spiral shape will emerge by analogy to the steady-state growth shapes seen for bulk crystal shapes (in both two and three dimensions).

The simplest expression for a spiral shape evolution model is based on following the length of each edge in time,(2.12)where *M* is the number of spiral edges that could exist in the first turn and *l*_{i} and are the length and tangential velocity of edge *i*, respectively. The tangential velocity is determined by equations (2.6)–(2.8). This spiral shape evolution model tracks the shape of one turn of a spiral from its initial formation as it evolves outwards away from the dislocation. Thus, in order to generate the full spiral from the set of *l*_{i}, equation (2.12) must be solved for many turns of the spiral and they must be placed end to end with the corresponding angle between each of the spiral sides. Nonetheless, it is only necessary to follow the dynamics of one turn (i.e. one revolution) of edges as they evolve away from the dislocation because each additional turn will follow identical dynamics while lagging in time by an integer multiple of the rotation time, *τ*. This spiral evolution model has no apparent steady state; however, by non-dimensionalizing the variables, steady-state features of the model are revealed. Gadewar & Doherty (2004) derived a similar set of equations for the two-dimensional crystals growing in time. The characteristic quantities used to define dimensionless variables for our system need to be chosen appropriately for a growing spiral. The characteristic length is chosen to be the perimeter of the turn of the spiral being considered, (the total length of all of the edges in one set *M*). In order to monitor the spiral evolution, a set of edges *M* is tracked as they move away from the dislocation. As the spiral moves away from the dislocation, *L* continually increases. The characteristic velocity is defined as the sum of all of the tangential velocities of all of the edges in that turn, . Thus, the dimensionless length is defined by *x*_{i}=*l*_{i}/*L* and the dimensionless velocity is defined by .

Using these dimensionless variables, the dynamic model is reformulated as follows:(2.13)and(2.14)A dimensionless warped time can then can be defined, , and equation (2.13) can be rewritten in fully dimensionless form,(2.15)The steady state of this system corresponds to the condition when each of the states of the model, *x*_{i}, is at a value that does not change with time. The steady state occurs when(2.16)When the step velocities (*v*_{i}) are constant in time, then , and hence *u*_{i} are also constant beyond stage (ii) of each edge's growth. Under these constant relative velocity conditions, the model is a linear system of ordinary differential equations, with all the eigenvalues equal to −1. It follows that the steady-state spiral shape (which corresponds to the shape of the spiral far away from the dislocation) is unique and stable.

The steady-state nature of the spiral shape is demonstrated for the spiral shown in figure 2*a*. Figure 3 shows the dimensionless length, *x*_{i}, of each of the three characteristic spiral edges (1–3), from their initial appearance at the dislocation and as they move outward away from the spiral. In the initial stages of the first rotation, the tangential velocities transition through their three stages of growth; thus, the curves are not smooth throughout that process. Far away from the dislocation, the dimensionless length approaches the steady state (*u*_{1,ss}=*x*_{1,ss}≃0.027; *u*_{2,ss}=*x*_{2,ss}≃0.188; *u*_{3,ss}=*x*_{3,ss}≃0.285). This highlights the stable steady state that all spirals characteristically evolve towards as they emanate from their dislocations. This model for spiral shape is general, in that it can also be used for more complex velocity profiles to determine the growth rate of the crystal face.

### (c) Spiral rotation direction

For an organic molecular system, the exposed chemical moieties along each edge of a given face are often different from one another. Thus, the edges have different velocities and critical lengths. Since the spiral direction will dictate which edge velocity is paired with the corresponding next edge's critical length, the growth rate resulting from the clockwise and anticlockwise spiral rotation directions can be different (see equations (2.3) and (2.5)). Additionally, different rotation directions could result in the disappearance of different edges. While there could be some preference to the direction of the spiral based upon the relative number of dislocations in each of the rotation directions, the actual growth rate of the face only depends on the most active spiral direction (Frank 1949). This is because the most active, or fastest growing, spiral direction will overtake the slower growth spiral direction. Thus, the characteristic spiral time of each spiral direction is calculated and the smaller *τ*, corresponding to the more active spiral (faster growth rate), will determine the normal growth rate of that face.

Figure 4 shows the spiral shape for each of the two rotation directions for a growth spiral. For this case, consider three independent edges (*h*, *j* and *k*) and thus six possible spiral sides. They have velocities and critical lengths (in arbitrary units) of *v*_{h}=2.0, *v*_{k}=1.0 and *v*_{j}=1.5, and *l*_{c,h}=2.0, *l*_{c,k}=3.0 and *l*_{c,j}=2.5, as well as angles *α*_{h,k}=*α*_{k,h}=45°, *α*_{k,j}=*α*_{j,k}=65° and *α*_{j,h}=*α*_{h,j}=70°. In both the clockwise and anticlockwise cases, we begin with side 1 as *h*. Thus, the clockwise spiral of figure 4*a* has its sides 1 and 4 as *h*, 2 and 5 as *k*, and 3 and 6 as *j*, giving velocities *v*_{1}=*v*_{4}=2.0, *v*_{2}=*v*_{5}=1.0 and *v*_{3}=*v*_{6}=1.5, with critical lengths *l*_{c,1}=*l*_{c,4}=2.0, *l*_{c,2}=*l*_{c,5}=3.0 and *l*_{c,3}=*l*_{c,6}=2.5 and angles *α*_{1,2}=*α*_{4,5}=45°, *α*_{2,3}=*α*_{5,6}=65° and *α*_{3,4}=*α*_{6,1}=70°. The anticlockwise spiral of figure 4*b* has its sides 1 and 4 as *h*, 2 and 5 as *j* and 3 and 6 as *k*, giving velocities *v*_{1}=*v*_{4}=2.0, *v*_{2}=*v*_{5}=1.5 and *v*_{3}=*v*_{6}=1.0, with critical lengths *l*_{c,1}=*l*_{c,4}=2.0, *l*_{c,2}=*l*_{c,5}=2.5 and *l*_{c,3}=*l*_{c,6}=3.0 and angles *α*_{1,2}=*α*_{4,5}=70°, *α*_{2,3}=*α*_{5,6}=65° and *α*_{3,4}=*α*_{6,1}=45°. In this case, the dominant growth spiral (thus the one that determines the face's normal growth rate) is the anticlockwise spiral in figure 4*b*, since its rotation time (*τ*_{anticlock}≃8.80) is less than that of the clockwise spiral in figure 4*a* (*τ*_{clock}≃9.16).

## 3. Predicting faces, edges, edge height, critical length and step velocity

In order to implement equations (2.3) and (2.5) to determine the crystal growth rates, we now need to predict which faces will be flat (F) on the crystal, as well as the edge directions and their critical lengths, step velocities and heights for each of the F crystal faces.

### (a) Face and edge determination

The number of strong periodic bond chains (PBCs) present on a particular crystal face helps to determine the mechanism by which that crystal face will grow. On faces with zero PBCs (kinked or K faces), or one PBC (stepped or S faces), the face tends to grow by the incorporation of molecules throughout the crystal surface. On faces where there are multiple PBCs, so termed F faces, growth occurs via a layer-by-layer mechanism where molecules incorporate along specific edges on the crystal surface, which are present owing to (i) spirals emanating from screw dislocations terminating at the surface or (ii) the formation of the two-dimensional nuclei on the surface. Since molecules incorporate throughout the surface on S and K faces, the growth rates of these faces are fast relative to the growth of F faces and are often limited by transport processes. Thus, we assume that all of the S and K faces grow much (10–1000 times) faster than the F faces (Lovette *et al*. 2008). Normally, the steady-state growth shape of a crystal will be bounded by slow-growing F faces. However, in some cases, particular directions may have no F faces to bound that direction and needle- or plate-like crystals will result.

At low levels of supersaturation, the spiral growth mechanism is dominant (Ohara & Reid 1973). In order to implement a spiral growth prediction, it is necessary to determine the edges that form the sides of the spiral. We assume that each of the PBCs forms a potential edge of the spiral and is included in the set *M*. These edges are relatively slow growing since the number of kink sites along them is far fewer than along any of the other directions on the face. In other words, edges running in any other direction on the face are expected to move with a velocity faster than their inclusion velocity (*v*^{incl}) and do not need to be included in the growth rate prediction.

### (b) Spiral edge height

The spiral edge height is determined by the distance between layers of growth. This corresponds to the smallest magnitude of the dislocation's Burgers vector that can provide for layerwise growth of the crystal. While some larger dislocations have been shown to exist on some crystal faces (Vekilov *et al*. 1992), each face is expected to be dominated by the smallest size Burgers vectors since the energy of the dislocation increases quadratically with the magnitude of the Burgers vector (*E*∼|** b**|

^{2}). Additionally, if larger sized dislocations are present (say twice the height of the smallest Burgers vector) when the steps flowing from them intersect with those from the smallest Burgers vector, those steps of twice the height will be reduced to half the size of the steps of the smallest height.

The height corresponding to the smallest magnitude Burgers vector is determined by first selecting a slice of the plane with a thickness of one interplanar spacing and determining whether it contains two or more PBCs. It is possible that the actual height of the step will be a fraction of the interplanar spacing. Therefore, the slice is then divided into half (or thirds depending on the symmetry). If each of the new slices contains growth units and is an F plane, then the division is performed again. This division continues until the newly divided slices do not each contain a set of growth units that form an F plane. The step height is then given by the thickness of the slice that corresponds to the last division, where each of the divided slices contained growth units that formed an F plane. In some complex systems, consecutive growth layers of a crystal face can have alternating properties, including rotation time. In these cases, it is necessary to evaluate the number of turns corresponding to one full grouping of the growth pattern. Additionally, as mentioned previously, Burgers vectors of magnitude greater than one have been seen. These situations can lead to potential step splitting or step doubling (Land 1997). The explicit modelling details of systems with these properties are not considered here and will be the subject of a future paper.

### (c) Spiral edge critical length

Predicting the critical length of a spiral edge requires the determination of the length below which the edge has a velocity of zero and above which the edge moves at a positive normal velocity (see equation (2.4)). While the growth of a spiral edge is a kinetic process, the determination of this critical length is often treated with thermodynamic theories. Kinetic theories for the determination of a critical length exist (Voronkov 1970); however, the physical properties necessary to implement them are not easily accessible or quantitatively predictable. On the other hand, the physical properties required for the thermodynamic models are more readily and rapidly attainable. Thus, thermodynamic methods are used here to estimate critical lengths.

One method, which has been used by Chernov (1984) and Winn & Doherty (1998) to determine the critical length of a spiral edge, is to predict the size of a two-dimensional critical nucleus using a classical Gibbs–Thomson approach. The critical length of a spiral edge is then assumed to be equal to the critical length for that same edge on an equivalent hypothetical two-dimensional nucleus. For this calculation, the critical length is determined by performing a free energy balance on the system and asking whether adding another molecule to that nucleus would increase or decrease the free energy of the system.

Rather than assuming that the critical length of a two-dimensional nucleus corresponds to the critical length on a spiral edge, we propose to determine the critical length using a free energy balance on a single spiral edge. This is similar to the approach proposed by DeYoreo and co-workers (Thomas *et al*. 2004). The critical length of edge *i* is then determined by asking whether adding a new row of molecules to this edge at its current instantaneous length either increases or decreases the free energy of the system (figure 5). If adding a row of molecules of length *l*_{i} increases the free energy of the system (Δ*G*>0), the edge is assumed to be unable to grow, while if adding a row of molecules of length *l*_{i} decreases the free energy of the system (Δ*G*<0), the edge is assumed to be able to grow. Thus, the critical length is the length at which the change in free energy caused by adding a new row of molecules to the system is zero.

Consider the case where edge *i*−1 is growing and continually lengthening edge *i*. To determine the critical length of edge *i*, consider the change in total Gibbs free energy (Δ*G*) when a new row of molecules is added to edge *i* of the spiral (figure 5), given by(3.1)where *n* is the number of molecules that were added in the new row; Δ*μ* (dimensions of energy per molecule) is the difference between the chemical potential of the molecules in solution and in the crystal; *γ*_{j} (dimensions of energy per area) is the specific surface energy of each new area being exposed; and *A*_{j} is the size of each new area being exposed. We assume that the surface energy along the length of edge *i* after the new row *i* is added is energetically identical to the surface energy along edge *i* before the new row is added. For this to be true for complex molecular systems (e.g. non-centrosymmetric molecules), multiple rows of molecules may need to be added. Thus, cancels from the calculation, and equation (3.1) can be rewritten as(3.2)where *h* is the height of the spiral edge; *a*_{p,i} is the distance the edge propagates with the addition of a row on edge *i*; *l*_{i} is the length of edge *i*; *V*_{m}/*N*_{A} is the molecular volume; is the edge energy (energy per area) of edge *i*−1; and *a*_{e,i} is the distance between molecules along edge *i* (figure 5). This can then be simplified to(3.3)where is the edge energy (energy per molecule) of edge *i* (; ). A characteristic plot for the free energy change for the addition of a new row is shown in figure 6. The critical length where Δ*G*=0 is given by(3.4)Note that the critical length of edge *i* depends on the edge energy of the adjacent edges of the spiral (*i*−1 and *i*+1). The distance between each molecule on edge *i* (*a*_{e,i}) is a simple function of the geometry. The chemical potential driving force (Δ*μ*) is assumed to be isotropic, thus is equal on each edge of each face (variables that are isotropic have the same value on each edge of each face and do not affect the relative rates necessary for shape prediction since they cancel as common factors). Critical lengths are expected to be in the range from 10 to 1000 nm. The only terms that remain to be predicted to determine the critical lengths are the edge energies per molecule, , for each edge of the spiral.

### (d) Spiral edge velocity

A constant competition exists between the attachment of molecules to the crystal and the detachment of molecules from the crystal. The most important area of molecular incorporation is along step edges. When the rate of attachment is greater than the rate of detachment, a crystal edge (step) will grow; however, when the rate of detachment is greater than the rate of attachment, the step will dissolve. When the rate of attachment is equal to the rate of detachment, the edge is in equilibrium. This attachment or detachment takes place as a four-step process. For attachment, the solute first diffuses from the bulk solution through a boundary layer to the crystal surface. Next, the solute molecule adsorbs onto the surface, then diffuses across the surface and finally is incorporated into the lattice. The detachment process is identical, but occurs in reverse. For organic molecular crystals grown in solution, the molecules often diffuse from the bulk directly to the locations on the surface where they incorporate (i.e. surface diffusion is a fast process). Additionally, for solution growth, the growth rate is often limited by the incorporation of molecules into the lattice (Chernov 1984). Thus, we use a model which assumes that lattice incorporation is the rate-limiting step. The most important site for the incorporation of molecules into the crystal is the kink site. In this paper, we focus on the growth velocity of crystals where a kink site can be defined classically as a ‘half-crystal position’ (Kossel 1927; Stranski 1928). The half-crystal position is defined as a location where exactly half the bonds available to that site are exposed to the solution, and the other half are interacting with the solid. This can only be valid for the case where the bond network surrounding each molecule is centrosymmetric. Determination of the step velocity for more complex bond network systems will be the subject of a future paper.

An expression for the step velocity of the *i*th edge where classical kinks along step edges are the primary sites of incorporation and disincorporation from the crystal, and where kink integration is rate limiting is given by(3.5)where *a*_{p,i} is the distance of propagation (figure 5) of edge *i* and *ν* is the frequency of attachment and detachment attempts (Markov 2003). The average spacing between kinks is *δ* and the molecular spacing along the step is *a*_{e}; thus, the quantity is the probability of finding a kink at any location along the step edge. The concentration of the solute in the solution is *C*, the saturation concentration is *C*_{sat} and *V*_{m} 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 (*ν*), the difference in concentration between the bulk and saturation (*C*−*C*_{sat}) and the exponential terms are each assumed to be isotropic, and thus all cancel in the calculation of relative growth rates and therefore have no effect on crystal shape. (For complex systems, the half-crystal position does not exist; therefore, these isotropic assumptions are not necessarily all valid. For such systems, the step velocity calculation may require a more specific face-dependent determination of these variables.) The kink density, or the probability of finding a kink site, was first derived by Burton *et al*. (1951). Similarly, we assume that the rearrangement of kinks occurs on a time scale that is faster than that of kink integration; therefore, the kinks are assumed to be in their most probable distribution and are Boltzmann distributed. Thus, the expression for the kink density when positive and negative kinks have different kink energies can be expressed as(3.6)where and are the energies required to form a positive and negative kink, respectively (energy per molecule). For the case of equal positive and negative kink energies , this reduces to the familiar expression(3.7)The final form of the step velocity is then given by(3.8)where *a*_{e}/*δ* is given by equation (3.6). Since *a*_{p,i} is readily calculated from the crystal geometry (i.e. crystallography), it only remains to predict the kink energy to determine the relative step velocities.

## 4. Interfacial energies

For crystal growth from solution, the edge and kink energies used in determining critical lengths and step velocities should be the interfacial energies between the solid and solvent. Using the above methods, these interfacial energies are the key variables remaining to be predicted for shape calculations. For this work, we use a classical method of interfacing two surfaces by using known and calculated data on the pure solute and solvent to determine the interfacial energies. In order to account for interfacial interactions in this way, we first estimate the solid and solvent interactions independently. Then, we use these values to estimate the interfacial energy between them.

To calculate the solid component of the edge energy, each of the solid broken bond energies that are exposed at the edge for a single molecule is summed component by component. In other words, the dispersive components to each of the broken bond energies exposed are summed together resulting in the solid-phase dispersive edge energy , and the Coulombic components to each of the broken bond energies are summed together resulting in the solute side Coulombic edge energy . The solute phase component to the kink energy has previously been shown to correspond to the broken bond energy of the intermolecular interaction of the molecules parallel to the edge (Winn & Doherty 1998). Thus, the dispersive and Coulombic components to the solute phase kink energy are simply the dispersive and Coulombic broken bond energies of that interaction. These quantities are each in the units of energy/molecule; however, the interfacial calculation requires the quantities to be on a per area basis. Thus, each of the above components is then recalculated on a per area basis (i.e. ).

In order to estimate the solvent side surface energy, we use a technique based upon solubility parameters, which has been proposed by Kaelble (1971) and is given by(4.1)where is the molar volume of the solvent; *δ*_{z} is the solubility parameter associated with interaction component *z* (*z*=Coulombic, h-bond donating, h-bond accepting, dispersive, etc.); and *f* is a fractional factor of the solvent internal energy displayed at an interface. Data for these variables for a variety of common solvents, in particular for the molar volume and solubility parameters, are available from a variety of sources (Kaelble 1971; Barton 1975). Sometimes, only a single solubility parameter for hydrogen bonding is reported and this then needs to be distributed into both accepting and donating components. This apportionment is carried out using , where the amount of donating and accepting fractions are estimated based upon the molecule (i.e. for equally donating and accepting solvents such as water or ethanol ).

Finally, the solute and solvent side components are used to determine an estimate for the energy required to form an interface between the solute and the solvent. This is done using a classical method where the interfacial energy between two surfaces can be estimated as the sum of the cohesive energy of each surface (the energy required to create the interface of each of the surfaces in a vacuum) minus the work of adhesion (the energy relieved by putting the two surfaces together). In general, the interfacial energy between two surfaces (a,b) is thus given by(4.2)where *γ*_{a} and *γ*_{b} are the interfacial energies of surfaces a and b (the cohesive energy) and *W*_{ad} is the work of adhesion. Furthermore, the work of adhesion can be estimated using the geometric mean rule (Israelachvili 1992). This method has been extended by Winn & Doherty (2002) to also focus on hydrogen bonding.

Here, we will use a generalized version of these equations that includes all of the existing interactions in both the solute and the solvent. These interactions, in general, consist of both a dispersive component (due primarily to van der Waals interactions) and a Coulombic component (which can stem from a variety of phenomenon such as an ionically charged species, hydrogen bonding components or electrostatic interactions resulting from pi bonding, particularly in aromatic compounds). Since the cohesive energy terms are those of the pure component, all (dispersive and all Coulombic) of the interactions from both the solute and the solvent need to be included. However, the work of adhesion only should include the appropriate terms where the solute and the solvent can interact.

For example, for systems where either the solute or the solvent has no Coulombic nature at all, the work of adhesion only results from the dispersive terms and we use (note that either or both of the Coulombic components of the cohesive energy in this expression will be zero)(4.3)Alternatively, for a system where significant hydrogen bonding is present, we use the following:(4.4)where the donating and accepting components to the hydrogen bonding character are distributed from the total Coulombic energy using the method described for those solvents in the previous paragraph.

In other systems where substantial pi bonding and aromatic rings are present, we use the following:(4.5)In this case, we have included a term in the work of adhesion to account for the negatively charged pi bonds interacting with the positively charged hydrogen bond-donating groups of the solvent. (If a solvent such as benzene were used and the solute had hydrogen bond-donating groups, the solute and the solvent in the final term would be reversed.)

In each case, the first four terms represent the work of cohesion and the remaining terms represent the work of adhesion from each of the components at the interface. In some cases, for more complex molecular crystals, the Coulombic character may stem from a combination of components (e.g. a combination of pi bonds and hydrogen bonding). For these cases, a more detailed account for each of these interactions may be necessary.

### (a) Crystal shape determination

Using the calculated relative growth rates outlined above, the steady-state crystal shape is given by the Frank–Chernov condition (Frank 1958; Chernov 1963) as(4.6)where *R*_{i} is the relative growth rate of face *i* and *x*_{i} is the relative perpendicular distance from the centre of the crystal to face *i*. Additionally, faceted crystal shape dynamics in both growth and dissolution can be estimated using faceted shape evolution methodologies (Zhang *et al*. 2006; Snyder & Doherty 2007).

## 5. Results

### (a) Naphthalene

The first example system used to demonstrate our crystal shape prediction methodology is naphthalene. The morphology of naphthalene has been studied previously using both the attachment energy method (Grimbergen *et al*. 1998) and the molecular simulation (Cuppen *et al*. 2004). Each of these methods was accurate at reproducing the shape of naphthalene grown from the vapour; however, they did not account for solvent effects. Experimental results for the shape of naphthalene grown from multiple solvents (cyclohexane and ethanol) are also available (Grimbergen *et al*. 1998). Thus, this system is a good choice for demonstrating the predictive power of our model, including its ability to predict solvent effects.

Naphthalene crystallizes in the monoclinic space group P21/a with the following lattice parameters: *a*=8.266 Å; *b*=5.968 Å; *c*=8.669 Å; *α*=*γ*=90°; and *β*=122.92 (Pawley & Yeats 1969). There are two molecules per unit cell and they are located at the unit cell positions (0,0,0) and . The energetic interactions between the molecules in the solid are calculated using HABIT (Clydesdale *et al*. 1996), with the force field of Scheraga (Némethy *et al*. 1983) assuming a monomer growth unit. The atomic charges for the HABIT calculation are obtained using Gaussian (Frisch *et al*. 2004). The strong intermolecular interactions in the solid state (those with an energy greater than *kT*) are given in table 1. Owing to the crystal and molecular symmetry, each of these interactions repeats to form one of the six strong bond chains that are used to determine the growth rate of the crystal faces. (A bond chain with an asterisk identifier is related to the bond chain of the same letter by symmetry.) Based upon the interactions listed in table 1, the families of the F faces of naphthalene are {001}, , {110}, and {200}. Since experimentally obtained crystal shapes are available for naphthalene grown from ethanol and cyclohexane, these two solvents are chosen for comparison with the predicted morphologies. While no direct evidence exists for the spiral mechanism on naphthalene, atomic force microscopy (AFM) results have shown it to be the dominant mechanism for some cases of growth of a similar molecule, anthracene (Cuppen *et al*. 2004). Thus, for each of the F faces, we implement the spiral growth rate calculation. The solvent data required to implement the shape prediction model is provided in table 2. For the case of growth in ethanol, the solute–solvent interactions are determined using equation (4.5) because of the polar ethanol molecule's interactions with the directional pi bonding of naphthalene.

Table 3 lists each of the F faces of naphthalene, the PBCs present on that face and the predicted relative growth rates of each face grown in both ethanol and cyclohexane. Detailed results for all the intermediate data (kink energies, edge energies, kink probabilities, etc.) for growth in cyclohexane and ethanol are given in the electronic supplementary material. Figure 7 shows the progression of the calculation. We begin with a known crystal structure (figure 7*a*). Then the intermolecular interactions in the solid state are calculated, which give rise to the PBCs in the system (figure 7*b*). These bond chains then correspond to the edges used for the predicted spiral shapes (figure 7*c*). Finally, the growth rates corresponding to these spirals lead to a crystal shape (figure 7*d*). Note that the spiral shapes indeed evolve towards a steady-state growth shape as they move away from the spiral centre. The crystal shapes based upon the growth rates in table 3 for growth in ethanol and cyclohexane are shown in figure 8, alongside the redrawings of the experimental results from Grimbergen *et al*. (1998). The predicted results using our spiral growth model are in very good agreement with the experimental results.

### (b) Anthracene

The morphology of anthracene has also been studied previously using both the attachment energy method (Docherty & Roberts 1988; Grimbergen *et al*. 1998) and the molecular simulation (Cuppen *et al*. 2004). These methodologies were accurate at reproducing the shape of anthracene grown from the vapour; however, they did not account for solvent effects. Experimental results have been previously reported for the growth of anthracene in an unspecified alcohol (Groth 1919). However, owing to unreported details associated with these results, we have performed our own experiments by slow evaporation of anthracene from a 2-propanol solution. The crystals were grown by preparing a slightly undersaturated solution of 2-propanol at room temperature and adding approximately 20 ml of it to the bottom of several small Petri dishes. The dishes were then each partially covered to varying extents and were let stand until all of the solvent was evaporated. The crystals grown in the mostly covered Petri dishes were well formed, and the results from these are used for comparison with growth predictions.

Anthracene crystallizes in the monoclinic space group P21/a with the following lattice parameters: *a*=8.562 Å; *b*=6.038 Å; *c*=11.184 Å; *α*=*γ*=90°; and *β*=124.70 (Mason 1964). There are two molecules per unit cell and they are located at the unit cell positions (0,0,0) and . The interactions between the molecules in the solid were predicted using HABIT (Clydesdale *et al*. 1996), with the force field of Scheraga (Némethy *et al*. 1983) assuming a monomer growth unit. The atomic charges for the HABIT calculation were obtained using Gaussian (Frisch *et al*. 2004). The strong interactions in the solid state (greater than *kT*) are given in table 4. Again, each of these interactions repeats to form one of the six strong bond chains that are used to determine the edges on each of the F faces and hence the growth rate of the crystal faces. The spiral growth mechanism is again used for growth rate determination of the F faces, and thus is modelled using the methodology described earlier in this paper. AFM results for vapour-grown anthracene crystals have been previously reported (Grimbergen *et al*. 1998), which showed spiral growth.

Based upon the interactions listed in table 4, the families of F faces of anthracene are {001}, , {110}, and {200}. The data on 2-propanol required to implement the shape prediction model is provided in table 2. Table 5 lists each of the F faces of anthracene, the molecular interactions corresponding to each bond chain present on that face and the relative growth rates of that face grown from 2-propanol. As was the case for naphthalene, the solute–solvent interactions are estimated using equation (4.5) owing to the pi-bonded rings in the compound. Detailed results for all of the intermediate data (kink energies, edge energies, kink probabilities, etc.) can be found in the electronic supplementary material.

The predicted crystal shape based upon the growth rates in table 5 for growth from 2-propanol is shown in figure 9, alongside several of the experimentally grown crystals. The experimental results show some degree of variation owing to the varied conditions of growth during solvent evaporation. Nonetheless, the predicted results using our spiral growth model are in good agreement with the typical shape of the experimentally grown crystals. Moreover, our experimental results are also similar to the shape reported by growth from alcohol (Groth 1919).

### (c) α-Glycine

Finally, we report the predicted results for α-glycine grown from water. The morphology of α-glycine has been studied previously using kinetic crystal growth methods (Bisker-Leib & Doherty 2003), but here we report successful morphology predictions without the need for expert decisions during implementation of the model. Experimental crystal morphologies of α-glycine grown in water are available from multiple literature sources (Boek *et al*. 1991; Poornachary *et al*. 2007).

α-Glycine crystallizes in the monoclinic space group P21/n with the following lattice parameters: *a*=5.1054 Å; *b*=11.9688 Å; *c*=5.4645 Å; *α*=*γ*=90°; and *β*=111.697 (Jönsson & Kvick 1972). There are four molecules per unit cell, and the growth unit is taken to be a hydrogen bonded cyclic dimer incorporating the zwitterionic interactions. The dimer growth unit can be identified from the solid-state interactions between the molecules, and has also been confirmed both by diffusion (Chang & Myerson 1986) and AFM (Carter *et al*. 1994). The two growth units are located at the unit cell positions (0,0,0) and . The interactions between the growth units in the solid were predicted using HABIT (Clydesdale *et al*. 1996), with the force field of Scheraga (Némethy *et al*. 1983) based upon the dimer growth unit. The atomic charges for the HABIT calculation were obtained using Gaussian (Frisch *et al*. 2004). The strong interactions in the solid are given in table 6. They repeat as four bond chains, two of which have identical interaction energies.

Based upon the strong interactions listed in table 6, the families of the F faces of α-glycine are {020}, {110} and {011}. α-Glycine is commonly grown from water and the solvent data required to implement our model is given in table 2. Since the growth unit consists of the hydrogen bonded cyclic dimer that includes the zwitterionic effects, all of the hydrogen bonding and zwitterionic capabilities of the solute are assumed to be satisfied internally to the growth unit (*γ*_{solu,acc}=0; *γ*_{solu,don}=0). The solute–solvent interactions are calculated using equation (4.4) (which reduces to equation (4.3) for this case of *γ*_{solu,acc}=0; *γ*_{solu,don}=0). Table 7 lists each of the flat faces of α-glycine, the molecular interactions corresponding to each bond chain present on that face and the relative growth rates of that face based on the spiral growth mechanism for growth in water. The detailed results for all of the intermediate data (kink energies, edge energies, kink probabilities, etc.) can be found in the electronic supplementary material. The crystal shapes based upon the growth rates in table 7 for growth from water are shown in figure 10*a*, and are very similar to the experimental results from Poornachary *et al*. (2007), shown in figure 10*b*.

## 6. Conclusions

We have reported on the next generation of kinetically based crystal growth models for the spiral growth mechanism. The evolution of spiral shape and the determination of which edges to include in the spiral shape calculation have been mathematically determined. This methodology has now advanced to become completely automated using a Matlab program, which can take crystallographic, solvent and molecular interaction information as input data and use the spiral growth methodology to systematically predict the steady-state crystal shape of organic molecular solids grown from solution.

While the model has been automated, and it appropriately accounts for crystallographic orientation, solvent and spiral dynamics, it also has some remaining limitations that are areas of future research. The determination of edge velocity has been derived only for systems where classical kinks (half-crystal positions) are present. Many crystal systems do fall into this category; however, many complex molecular organic crystals such as pharmaceutical molecules do not. Thus, a methodology for a more complete description of the edge velocity for systems without classical kinks is necessary. Additionally, for those more sophisticated systems, a more complex methodology for edge and face determination will also be necessary. For this work, the resulting relative crystal growth rates that determine the steady-state crystal shape are currently independent of supersaturation. However, it is well known that some crystals have supersaturation-dependent relative growth rates when grown from solution (e.g. paracetamol grown from water; Shekunov & Grant (1997)). Methods to account for supersaturation dependencies on relative growth rate for the spiral growth model, and methods to appropriately determine the transition from spiral growth to the two-dimensional nucleation and growth are also desired.

One hallmark of our approach is that it is modularized in nature. When higher fidelity methods to determine edge velocities, critical lengths, solute–solvent interactions or any other estimated values in the model are developed, they can be readily incorporated into the calculation procedure. Furthermore, the model can even be implemented when some of these values are taken from the experiments. The key idea is to base the crystal growth rates on the simplest physically realistic representation of the underlying molecular mechanisms, thereby linking crystal engineering with crystal chemistry. Using these strategies, we have successfully developed and implemented a rapid method for the prediction of crystal shape that serves as a basis for future developments in crystal shape engineering, including mixed solvents and additives.

## Acknowledgments

We are grateful for the financial support provided by Merck & Co. and the National Science Foundation (grant no. CBET-0651711).

## Footnotes

- Received June 5, 2008.
- Accepted November 28, 2008.

- © 2009 The Royal Society