Predictive Modeling of Supersaturation-Dependent Crystal Shapes

Article pubs.acs.org/crystal

Predictive Modeling of Supersaturation-Dependent Crystal Shapes Michael A. Lovette and Michael F. Doherty* Department of Chemical Engineering, University of California, Santa Barbara, California 93106-5080, United States W Web-Enhanced Feature * S Supporting Information *

ABSTRACT: Several molecular organic crystals are known to exhibit significant changes to their steady-state growth shapes when grown in solutions with different supersaturations. In this article, we describe a first-principles model for the prediction of supersaturation-dependent growth shapes that is capable of capturing this phenomenon. The model begins with the prediction of the growth shape at an arbitrarily low supersaturation, for which all faces present on the steady-state shape grow by spirals. The dominant growth mechanisms (spiral or two-dimensional nucleation) are then determined over a specified range of supersaturations, based on the surface energies of edges present on the faces of the low-supersaturation (spiral growth) shape. These energies are estimated from the set of strong intermolecular interactions contained within each face and solubility parameter data. Steady-state growth shapes are then predicted as a function of supersaturation. A case-study performed on naphthalene grown in ethanol is presented in order to demonstrate the implementation of this model. The predictions from this case study are in remarkable agreement with the results of prior experiments performed for this system.

■

previously developed for spiral growth.6 This model for 2D nucleation determines the relative growth rates of faces from the intermolecular interactions contained within the crystal and the interactions occurring across solid−solution interfaces. Specifically, nuclei are treated as faceted objects bounded in directions determined by the PBCs contained within their respective faces. This is a departure from the Wulff construction, which has conventionally been used to determine the shapes of nuclei, as the Wulff construction allows for nuclei that are bound by edges that are not coincident with PBCs. This was done because the edge energy calculation for PBC directions is inexpensive and can be implemented in an automated manner, in contrast to the unrestricted Wulff construction that requires more expensive calculations for nearly all systems (the exception being Kossel crystals, for which closed form expressions have been determined for edge energies in arbitrary directions). The article concludes with a case study performed for naphthalene crystallized in ethanol that provides detailed methods for estimating the parameters used by our model. The steady-state growth shapes we predict using this model are in agreement with the shapes obtained experimentally.10

INTRODUCTION The recognized impact of particle shape on downstream processing and performance characteristics for solution crystallization products has led to the development of various methods for growth shape prediction.1−6 The models developed by Doherty et al.1,2,6 predict crystal growth shapes from first-principles based on the determination of key anisotropic parameters for each face and are rapidly implementable, a characteristic that makes them attractive for process design. However, these models have assumed that the growth of all F-faces (faces with two or more periodic bond chains (PBC)s as defined in ref 7) on a specific crystal proceeds exclusively by the same mechanism. This assumption restricts their application to limited ranges of supersaturation, for which all F-faces grow by either spirals (low supersaturation) or two-dimensional (2D) nucleation (high supersaturation). Therefore, these models cannot be used to estimate shapes over regions of moderate supersaturation for which spiral growth is the dominant growth mechanism on some F-faces, while 2D nucleation is dominant on others. Furthermore, existing models do not provide a systematic method to determine their applicable ranges of supersaturation. These limitations render such models incapable of capturing the significant changes to growth shapes, which occur over narrow ranges of supersaturation for a variety of materials.8−11 In this article, we present a first-principles mechanistic model for shape prediction that accounts for the effect of supersaturation on growth shapes and retains the desirable characteristic of rapid implementation. The model allows for different F-faces on the same crystal to grow by either a spiral or a 2D nucleation mechanism, with the determination of mechanism carried out systematically on a face-by-face basis at each supersaturation investigated. Additionally, a new model for growth by 2D nucleation was developed in order to model growth by this mechanism in a manner that is consistent with the universal approach © 2012 American Chemical Society

■

GROWTH MECHANISMS The three most common mechanisms for crystal growth from solution are (1) layered growth by spirals, (2) layered growth through the formation of 2D nuclei, or (3) rough growth. The growth of a crystal in the direction perpendicular to a face will proceed by whichever mechanism provides the fastest rate. This rate, which is hereafter referred to as the perpendicular growth rate of a face, depends on supersaturation and intermolecular Received: July 7, 2011 Revised: December 3, 2011 Published: January 20, 2012 656

dx.doi.org/10.1021/cg200855p | Cryst. Growth Des. 2012, 12, 656−669

Crystal Growth & Design

Article

and in most cases will not be present on a crystal’s lowsupersaturation steady-state growth shape. On faces containing two or more PBCs (referred to as F-faces),7 kink sites are located along edges which are spaced throughout the face. Edges are defined as the intersections between successive layers on a face. The growth of an F-face requires the creation and spreading (which occurs through the outward normal motion of edges) of these layers and is referred to as layer-by-layer or layered growth. Kink sites are the energetically and kinetically favored sites for the incorporation of molecules into an edge. Mechanisms for kink formation have been treated elsewhere.14−17 The average number of molecules along an edge between successive kink sites is given by

interactions (solid−solid and solid−solution) present at the face (among other factors). For solution grown crystals, these interactions are manifested as the surface energies of edges and kink sites. The dependence of growth mechanism on supersaturation for F-faces is shown in Figure 1. The progression of

⎛ ϕkink ⎞ 1 ni̅ = = exp⎜⎜ i ⎟⎟ + 1 ae, i 2 ⎝ kbT ⎠

Figure 1. The dependence of growth mechanism on supersaturation illustrated for an F-face. Adapted with permission from ref 5. Copyright 2009 American Chemical Society.

x0, i

(1)

for edge i; where x0,i is the average distance between kink sites, ae,i is the distance along the edge occupied by a growth unit, ϕikink is the work required to create a disturbance along an edge (dimensions of energy per site), kb is Boltzmann’s constant, and T is the temperature (for reference kbT ≈ 4.14 × 1014 erg at 300 K). Equation 1 assumes the edge is in its equilibrium (most probable) configuration, and there is only a single “type” of kink site found along the edge. The assumption of a single type of kink, which is applied hereafter in this article, restricts this analysis to crystals containing only centrosymmetric PBCs. A similar analysis for crystals with multiple types of kink sites has recently been developed by Kuvadia and Doherty.18 Equation 1 is valid whether disturbances along the edge are created by thermal fluctuations, or the average separation between kinks is maintained by a one-dimensional (1D) nucleation mechanism with Si < n̅i−1, where Si = (Ci/Csat) − 1, and Ci, Csat are the concentrations of solute molecules at the edge and in a saturated solution, respectively.15,17 This expression for kink density was formed in the limit of ϕikink ≫ kbT, where all disturbances to an edge result in the formation of a kink site. This condition is not met when ϕikink ≲ 3kbT; in this case eq 1 can be modified.19 For edges not coincident with PBC directions (e.g., for instance, the black edges in Figure 3b or Figure 4b) ϕkink ≪ kbT and the separation between successive kink kink sites will trend toward a minimum of n(ϕ → 0) = 2.17,19 ̅ These edges will move in their outward normal directions at an exceedingly fast rate when compared to those coincident with PBCs. As the strengths (adjusted to account for solid-solution interactions) of the PBCs contained within a face increase, there is a diminishing likelihood of observing edges that are not coincident with PBCs on a spiral beyond its first turn6 or on a larger than critically sized nucleus. When the set of PBCs contained within a face are all weak (ϕikink ≲ kbT for all edges), round or elliptical spirals and 2D nuclei containing several edges that are not coincident with PBCs are expected as, roughly speaking, n̅i ≈ 2−3 for all edges on the face. In layered growth, edges spreading in their outward normal directions result in the addition of a new layer of molecules to the crystal in the direction of the face every τ units of time. Therefore, the perpendicular growth rate of the hkl face is given by

growth mechanisms with increased supersaturation has been observed experimentally using atomic force microscopy by Land and De Yoreo.12 Recently, Sleutel et al.13 have experimentally investigated the microstructural surface dynamics accompanying growth by 2D nucleation and rough growth and present a remarkably clear view of these processes. The maintainable supersaturation in a conventional batch cooling crystallizer is limited by the system-specific metastable zone width. This limitation arises from spontaneous bulk nucleation (creation of new crystals), which occurs at supersaturations beyond the metastable zone width; a phenomena that rapidly depletes the supersaturation. Furthermore, the growth rate of a face can be limited by the kinetics of surface integration processes and/or the rates of solute and heat transport toward or away from the face.5 Figure 1 can be interpreted as follows: the growth rate of the face as a function of supersaturation is given by the solid curves; the red line indicates the bulk transport (BT) limit; the black curves represent the mechanism-dependent surface integration limits;5 the points σ2D,hkl and σBT,hkl indicate the onset supersaturation for growth by 2D nucleation and bulk transport limited growth, respectively; and the vertical green lines marked MZ1 and MZ2 are two possible positions for the metastable zone boundary. In Figure 1, the bulk transport limit is drawn such that σ2D,hkl occurs at a lower supersaturation than σBT,hkl. This is not a requirement and depends on several system specific parameters including: the surface energies of edges and kink sites, temperature, the viscosity of the solution, and mixing.5 Two cases are presented in this figure. In case “MZ1”, the metastable zone boundary occurs at a lower supersaturation than σ2D,hkl; therefore, in this case the hkl face will typically grow exclusively by a spiral mechanism. In the case “MZ2”, the boundary occurs at a higher supersaturation than σ2D,hkl, and the dominance of growth by a 2D nucleation mechanism is feasible.

■

LAYERED GROWTH In the surface integration limit (as defined in ref 5), the perpendicular growth rate of a face depends on the net rate of incorporation into kink sites and the average distance between them. On faces containing zero or one periodic bond chain (referred to as K- and S-faces, respectively),7 kink sites are located throughout the face. These faces grow rapidly (in comparison to F-faces) by a rough mechanism at any supersaturation

Ghkl = 657

⎛h⎞ ⎜ ⎟ ⎝ τ ⎠hkl

(2)



Article

given by VmC = VmCsat(S + 1) for a supersaturated solution. Replacing VmCsat with VmCsat(S + 1), in eqs 5 and 6 and simplifying the resulting expression, τrow is given by

where h is the height of a layer; which is typically equal to the interplanar spacing in that direction or a small integer fraction or multiple thereof. The coverage time of a face, τ, depends on the growth mechanism, edge velocities, critical lengths, and velocity profiles. These features are determined by the supersaturation and hydrodynamic/thermal conditions prevalent during growth. Due to their importance in growth processes and therefore shape prediction: edge velocities, critical lengths, and velocity profiles are described in the following sections. Edge Velocities. The motion of an edge its outward normal direction occurs through the net incorporation of molecules into kink sites. The velocity of this motion is given by 1 vi = a p, i τ− row, i

n̅ (8) Sk− for a supersaturated solution. Therefore, the velocity of an edge in its outward normal direction is given by τrow =

⎛ k− ⎞ v = a pS⎜ ⎟ ⎝ n̅ ⎠ (9) For the solution growth of small organic molecules, substituting typical order of magnitude estimates of v = 6 (1000) nm/s, n̅ ≈ 10, ap ≈ 1 nm and S ≈ 0.1 (based on refs 9 and 5) in eqs 9 and 8, result in k− = 6 (105) s−1, and τrow = 6 (0.001) s. Edges with conditions corresponding to these typical values advance in their outward normal directions at a rate that is on the order of 1000 rows/s. Equation 9 is often expressed using the stepkinetic coefficient β, defined such that

(3)

where ap,i is the distance edge i moves with the incorporation of a new row of molecules into the edge, and τrow,i is the average time required for this incorporation to occur. The incorporation of molecules into kink sites proceeds in a self-sustaining manner. Kinks are unoccupied sites; therefore, an incorporation event into a kink site “destroys” the kink site in question, which following this event becomes classified as an edge site. However, with each incorporation event a kink site is formed at a neighboring position.17 The net effect of this process can be described as the “motion” of a kink. Assuming an edge is maintained in its most-probable state (through thermal roughening or 1D nucleation processes occurring at a faster pace than incorporation events), its displacement through the incorporation of a new row of molecules occurs on average once each kink site has moved ni̅ times. The time for this process to occur is given by

n τrow = + ̅ − j −j +

k+ + k− n̅ For systems with k+ ≫ k−, β can be approximated by β = ap

k+ (12) n̅ Therefore, eq 12 that is often used to define β, is most accurate for systems with low solubility (see eq 7). Expressions for k+ or k− can be constructed using reaction rate theory. Since Burton, Cabrera and Frank,14 transition state theory (ballistic motion over an energy barrier) has been applied almost ubiquitously17,21 to determine expressions for k+ and k−, resulting in

(4)

where j and j are the fluxes into and out of kink sites, respectively. The subscript i, which references a quantity to a specific edge, is dropped in eq 4 and throughout this section. In a saturated solution, these fluxes are given by20

j = VmCsatk

+

(11)

β ≈ ap

−

+

(10)

v = VmCsatSβ Applying this definition of β to eq 9 results in

⎛ ΔE‡ ⎞ k T ⎛ Z‡ ⎞ + ⎟⎟ k TST = b ⎜⎜ ⎟⎟ exp⎜⎜− 2πℏ ⎝ Z0 ⎠ ⎝ kbT ⎠

(5)

and

(13)

and

j− = (1 − VmCsat)k− +

(6) − k TST =

−

where k and k are the rate constants for the attachment and detachment of molecules into and out of kink sites, which have dimensions of inverse time, and Vm is the molar volume of the solute. The quantity VmCsat is equal to the probability of encountering solute molecules in a saturated solution. Alternatively, VmCsat can be thought of as the number density of solute molecules in a saturated solution. In a saturated solution these fluxes balance, and edges are stagnant with v = 0 and τrow = ∞. From this balance (i.e., j+ = j−), the relationship

k− VmCsat = + k + k−

⎛ ⎛ ΔE‡ + ΔW kink ⎞⎞ kbT ⎛ Z‡ ⎞ ⎜⎜ ⎟⎟ exp⎜⎜− ⎜⎜ ⎟⎟⎟⎟ 2πℏ ⎝ Z0 ⎠ kbT ⎠⎠ ⎝ ⎝

(14)

‡

where ℏ is the reduced Planck’s constant, Z , Z0 are the partition functions of the transition and initial states, respectively, ΔE‡ is the energy of the transition state, and ΔWkink is the work required to remove a molecule from a kink site. The ratio of Z‡/Z0 and kbT/2πℏ have often been grouped and the edge velocity expressed as5,14,16,22

⎛ ΔE‡ ⎞ ⎛ ap ⎞ ⎟⎟ vTST = VmCsatS⎜ ⎟υ exp⎜⎜− ⎝ n̅ ⎠ ⎝ kbT ⎠

(7)

(15)

Equation 15 is arrived at by combining eqs 10, 12, and 13 and substituting υ ≡ ((Z‡)/(Z0))((kbT)/(2πℏ)). The edge velocity calculated using transition state theory has a m−1/2 dependence on the mass, m, of the growth unit that originates in the vibrational component of the partition functions in eq 14. In a recent series of papers by Vekilov et al.,17,23−26

emerges between the number density of solute molecules in a saturated solution and the rate constants for the attachment and detachment of molecules at kink sites. In the (strict) absence of bulk transport limitations, the probability of encountering a solute molecule at the edge is 658



Article

In the absence of kinetic limitations, the same critical lengths are expected for equivalent edges on a spiral or a 2D nucleus, with the nucleus’ shape corresponding to the minimum free energy for a given volume. This shape is derived in the section Two-Dimensional Nucleation Thermodynamics. Within this limiting case, the lengths of individual edges can be obtained by applying the Gibbs−Thomson law (GTL), which for four-sided spirals or 2D nuclei each with 2-fold rotational symmetry yields

edge velocities for the proteins ferritin and apoferritin, which have significantly different masses (MW = 450 000 and 750 000 g/mol for apoferritin and ferritin, respectively, where MW is the molecular weight), were shown to be equal when grown in aqueous solutions. Because of the lack of dependence on the mass of the solute molecule and the linear relationship found between ln(v) and T−1 (indicating the presence of an energy barrier), the kinetics of kink integration for these systems were proposed to follow Kramer's theory in the spatial diffusion limit (diffusive motion over an energy barrier).23,27 For this case, the rate constant is + kSD

=

⎛ ΔE‡ ⎞ ⎟⎟ exp⎜⎜− Λ ⎝ kbT ⎠

⎛+⎞ a p−1⎜ ⎟ ⎝ ⎠

lc, i

vSD

lc, i = 2 (17)

⎧ 0 for (li < lc, i) ⎪ ⎪ ⎛ vi = ⎨ exp(σlc, i/li) − 1 ⎞ ⎪ v∞ , i⎜1 − ⎟ for (li ≥ lc, i) ⎪ exp(σ) − 1 ⎠ ⎝ ⎩

(20)

(21)

where v∞,i denotes the steady-state velocity (v∞,i = vi(li = ∞)) of edge i. Removing the restriction of consistency with the GTL, Chernov et al. obtained31 an implicit expression for vi/v∞,i, given by

a p, i ni̅

ae, i ⎛ ϕikink ⎞ ⎜ ⎟ σ ⎜⎝ kbT ⎟⎠

When the conditions applied in eq 19 are true, eq 19 is recovered as a specific instance of eq 20. It is worth noting that eq 20 is consistent with the expectation that directions not coincident with PBCs (for which ϕikink ≪ kbT is expected) will not require the attainment of a critical length before moving in their outward normal directions. For the remainder of this article eq 20 is used to determine critical lengths for specific edges. Readers are directed to ref 29 for further discussion of critical lengths in cases where the GTL does not apply. Velocity Profiles. Teng et al. determined the velocity profile for edges moving in their outward normal directions in a manner consistent with the GTL as30

On the basis of a narrow grouping of activation barriers for the growth of a diverse set compounds in aqueous solutions (centered around ΔE‡ ≈ 11kbT at 300 K), Petsev et al. further asserted the prevalence of this limit for the solution crystallization of small molecules.24 Whether ballistic or diffusive motion over an energy barrier most accurately represents the kinetics of kink integration, the determination of this barrier along an appropriate reaction coordinate, and a corresponding vibrational partition function or diffusivity, are the minimum requirements for the prediction of edge velocities in the absolute sense (i.e., being scaled to “wall-clock” time). These computations remain expensive for the crystallization of small organic molecules from solution. However, assuming isotropic rate constants (i.e., ki+ = k+ and ki− = k− for all i), the relative edge velocities can be predicted from the remaining terms in eq 9. The designation as “relative” indicates velocities that are not scaled by the frequencies of crossing an energy barrier and therefore lack an appropriate time scale. Using this approach, the relative velocity of edge i is given by

vrel, i =

(19)

where vm is the molecular volume (vm ≡ Vm/NA where NA is Avogadro’s number), γj ≠edge i is the edge energy per unit area of the edges adjacent to edge i, and αi±1 is the interior angle between adjacent sides.5,28 Lovette and Doherty29 determined the critical length of an arbitrary edge obeying the GTL as

(16)

where + is the diffusivity of the solute molecule over the barrier and Λ is a length scale containing the radius of curvature of the 1D potential energy surface E around ΔE‡,24 resulting in

⎛ ΔE‡ ⎞ 1 ⎛+⎞ ⎟⎟ = VmCsatS ⎜ ⎟ exp⎜⎜− n̅ ⎝ Λ ⎠ ⎝ kbT ⎠

⎛ edge ⎞ vm ⎜ γ j ≠ i ⎟ 2 = sin(αi , i ± 1) σ ⎜ kbT ⎟ ⎝ ⎠

⎞ ⎛ ϕedge + ϕedge ⎞−1⎛ yi i−1 i+1 ⎟ ⎜ ⎟ ⎜ =⎜ ⎟ ⎜ 2 1/2 ⎟ lc, i k T − (1 y ) b ⎝ ⎠ ⎝ ⎠ i l* × (π + 2 arcsin(yi )) + i lc, i li

(18)

Equation 18 enables the prediction of relative edge velocities from the work required to form a kink (ϕikink) and the geometry of the crystal. Methods for predicting ap,i and ni̅ have been developed elsewhere.1,6,18,19 Critical Lengths. The critical length of an edge defines the average length required before it begins to move in its outward normal direction, which is to say that on average vi = 0 for li < lc,i, where li and lc,i are the actual and critical lengths of the edge, respectively. For edges on a spiral, this length is attained through the outward normal motion of the preceding edge in a process further described by Snyder and Doherty.6 In contrast, the edges of a 2D nucleus have attained a critical length once the nucleus as a whole has reached a shape-dependent critical size. Below this critical size the dissolution of the nucleus is a “favorable” process, and therefore, below this size edges are on average unlikely to move in their outward normal directions.

(22)

where yi ≡ (vi)/(v∞,i), ϕiedge is the edge energy on a per molecule basis (defined as ϕiedge ≡ γiedges, where s is the surface area exposed by a molecule on edge i) and li* is a correction to the critical length to account for the stochastic nature of the edge. This expression is bounded by 0 ≤ yi < 1. As ϕiedge ± 1 increases, the velocity profile given by eq 22 becomes a unit step function centered at li = lc,i, that is,

⎧ for (li < lc, i) ⎪ 0 vi ≈ ⎨ ⎪ ⎩ v∞ , i for (li ≥ lc, i) 659

(23)



Article

Layered Growth by Two-Dimensional Nucleation. For growth by 2D nucleation, τ2D denotes the coverage time, the time required to completely cover over a face with a single layer of nuclei. The perpendicular growth rate of a face growing by 2D nucleation is given by

This profile is referred to as the Voronkov velocity profile after ref 32. The different velocity profiles are shown in Figure 2. The region shaded for the GTL corresponds to 0.01 ≤ σ ≤ 1,

G2D =

approximating the typical range of operation for the industrial crystallization of small organic molecules. The GTL appears applicable in the range of kbT ≲ ϕedge i ± 1 ≲ 2kbT (for σ ≤ 1). To provide a frame of reference: for organic molecular crystals in solution, this region corresponds to edges typified by weak dispersive or favorable polar and hydrogen bonding interactions across the solid−solution interface. Such a system is exemplified by edges exposing hydrogen bonding groups in a strong hydrogen bonding solvent such as water. The Vornokov velocity profile (given by eq 23) becomes an increasingly more reasonable estimate as ϕedge i±1 ≳ 5kbT. For organic molecular crystals in solution, this region corresponds to edges typified by strong dispersive, polar, and/or hydrogen bonding interactions that are unfavorable in the given solvent. Such a system is exemplified by edges exposing strong hydrogen bonding moieties in a nonpolar solvent such as hexane. In this work, we apply the Voronkov velocity profile together with eq 20 to predict crystal shapes. Although for the cases investigated, there is a degree of inaccuracy introduced in using these expressions (in the absence of closed-form expressions for rotation times using a more applicable velocity profile) these inaccuracies are assumed to be isotropic. Layered Growth by Spirals. For spiral growth, τS is the characteristic rotation time, the time required for the innermost turn of the spiral to complete one revolution. The perpendicular growth rate of a face growing by spirals is given by

h τS

(24)

For N-sided polygonal spirals with edges adhering to the Voronkov velocity profile (eq 23), this time is given by N

τS =

∑ i=1

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

(26)

This time depends on nucleation rates and edge velocities. Conventionally, τ2D has been determined by assuming diskshaped nuclei with isotropic edge velocities, commonly referred to as the birth-and-spread model.22,33 However, this is an unlikely case for small organic molecules that crystallize almost entirely (∼95%) into lattices of orthorhombic or lower symmetry, that are bound by anisotropic van der Waals, hydrogen bonding, and electrostatic interactions.34,35 For anisotropic 2D nuclei, coverage times depend on properties defined by two separate shapes, namely, the “Wulff” and “kinetic” shapes. The Wulff shape is obtained from thermodynamic considerations and is used in determining the free energy barrier that must be crossed in the formation of a postcritical 2D nucleus. The kinetic shape is obtained from the relative edge velocities using the 2D analogue of the Frank-Chernov condition and determines the rate at which nuclei cover the face. The following events occur during the “life” of 2D nucleus (1) the agglomeration of molecules on the face leads to the formation of a critical nucleus typically in or near its Wulff shape; (2) the edges present on this nucleus that are not coincident to PBCs spread at an exceedingly fast rate in their outward normal directions becoming vertices between adjacent edges that are aligned in PBC directions; and (3) the spreading nucleus takes on its kinetic shape bounded by PBCs with its shape determined by the Frank-Chernov condition. Although there are several alternative models for crystal growth by 2D nucleation, each has its own shortcomings when used in shape predictions for this class of molecules. The model developed by Nielsen determined nucleation rates for several polygonal nucleus shapes by introducing a shape factor to account for the deviation from discs.36 While this provides increased accuracy for small organic molecules over models that assume disk-shaped nuclei, the model still assumes isotropic spreading rates. Lewis provided a framework for atomistic modeling of 2D nucleation; however, only the case of a Kossel crystal (simple cubic crystal with 6-fold symmetric intermolecular interactions) was considered, and the subsequent growth of nuclei (to postcritical sizes) was not addressed.37 The model developed by Winn and Doherty1 can be applied for any regular N-sided nucleus; however, it approximates the relative growth rates of different faces as Rrel,hkl ∝ 1/ΔGc,hkl where ΔGc,hkl is the free energy of formation of a critical nucleus on the hkl face. This model is applicable only for systems where all F-faces on the crystal grow by 2D nucleation and in the limit of nucleation (as opposed to spreading) rate limited growth kinetics. More recently, Cuppen et al.38 developed a model for growth by 2D nucleation with anisotropic spreading and shapedependent nucleation rates. However, their application of this model was limited to the specific case of rectangular nuclei. In this section, we generalize their approach for arbitrary polygonal shapes. Furthermore, growth by 2D nucleation has been simulated using Monte Carlo,39−41 kinetic Monte Carlo,42−46 and molecular dynamics.45,47−49 Because of the computational constraints,

Figure 2. Velocity profiles constructed using eqs 21, 22, and 23. The gray area corresponds to eq 21 for σ = (0.01 − 1). Solid curves edge / correspond to profiles obtained using eq 22 with the values of ϕi±1 kbT adjacent to each curve and assuming l ≈ 0.

GS =

h τ2D

(25)

where αi,i+1 is the angle between the vectors normal to sides i and i + 1 with αN+1,N = α1,N and N is the number of sides that remain on the spiral beyond its first turn.6 The scaling (i.e., functional dependence) of τs ∝ σ−2 and therefore GS ∝ σ2 for growth by spirals is evident from eq 25 coupled with the scaling of vi ∝ σ and lc,i ∝ σ−1 determined in the sections Edge Velocities and Critical Lengths for low supersaturations. 660



Article

been discussed with the following sections dedicated to the determination of J. Two-Dimensional Nucleation Kinetics. Two-dimensional nucleation is an activated process with an energy barrier given by ΔGc. Therefore, nucleation rates take the form

these simulations have been conducted: in the absence of solvent,39−41,43,44 for Kossel crystals,39,40,42,47 or at a single supersaturation.42,45,48,49 While such simulations have revealed important aspects of the physics of nucleation (often resulting in the improvement of mechanistic models), the computational expense incurred for organic molecules growing from explicit solvents limits their usefulness within a framework amenable to process design models.5 Following the derivation for rectangular-shaped nuclei provided by Cuppen et al.,38 a general expression can be developed for growth through the formation of N-sided convex nuclei with arbitrary edge velocities. Assuming edge velocities are constant, the area covered, ΔA, by a growing nucleus during the time interval Δt is given by

ΔA = 3f (vi , vit , N )(Δt )2

⎛ ΔGc ⎞ J = κ2D exp⎜− ⎟ = κ2Dξ ⎝ kbT ⎠

where κ2D is a kinetic prefactor and ξ is the dimensionless nucleation rate. Nielsen provides an approximation of κ2D ≈ + /h4, where + is the diffusion coefficient of the solute in solution.36 Making several “questionable” assumptions along the way (a characterization they themselves used), Ohara and Reid obtained

(27)

2v ̅ η12 ⎛ vm ⎞1/2 1/2 ⎜ ⎟ κ2D = σ π ⎝h⎠

f(vi,vit,N)

where relates the normal and tangential edge velocities, vi and vit, respectively, and the number of sides N, to the area covered by that nucleus per unit time and is given by

1 f (vi , vit , N ) = 6

∑ vivit (28)

The tangential velocity of an edge (i.e., the rate at which the edge lengthens) is determined from the velocities of adjacent edges in their outward normal directions, using the relation vit =

vi + 1 − vi cos(αi , i + 1) sin(αi , i + 1)

+

vi − 1 − vi cos(αi − 1, i) sin(αi − 1, i)

(29)

where vN+1 = v1 and αN,N+1 = αN,1. These expressions assume edges move with an “on/off” velocity profile consistent with eq 23. Edges with vit ≤ 0 are shrinking as a nucleus grows and will disappear during this process (becoming the vertices between other edges on the nucleus). These edges are not included in the expression for f. The condition vit = 0 results in the same critical velocity for the presence of an edge on a nucleus that was determined by Szurgot and Prywer.50 For regular polygons with vi = vj and αi = αj = 2π/N for i and j = 1, ..., N; f = 1/3cgvi2, where cg is a shape factor, such as those given by Nielsen.36 The generality of eq 28 allows for irregular shapes where differences in the critical lengths, angles between successive edges, and edge velocities are possible. The fractional coverage θ of a face growing by 2D nucleation at time t is given by38 6

θ = J f (vi , vit , N )t 3

ΔG = AP − VR

(31)

Substituting eq 31 into eq 2, the perpendicular growth rate of a face growing by 2D nucleation is given by

G2D = h(J f (vi , vit , N ))1/3

(35)

where the terms VR and AP are the free energy reward associated with the increase in volume of the more stable phase (“volume reward”) and penalty associated with the increase in the area of the interface between the two phases (“area penalty”), respectively. As defined, the quantities VR and AP are positive. For homogeneous 2D nucleation, this energy can be written in terms of the number of molecules in the nucleus, n, as

(30)

assuming all nuclei are the same size and shape, where J is the nucleation rate. The coverage time τ2D can be determined by solving for θ = 1, resulting in

τ2D = (J f (vi , vit , N ))−1/3

(34)

for κ2D, where v ̅ is the speed of a molecule striking a nucleus, and η1 is the number concentration of unbound molecules on the surface.22 Alternative expressions for κ2D are developed in ref 51. Using either approximation, Lovette and Doherty obtain κ2D = 6 (1025) cm−2 s−1, using the typical parameters for molecular organic crystals provided in Table 2 of ref 29. Moreover, the exponential dependence of nucleation rate on ΔGc make the accurate (in this article we argue that accurate can be ±10kbT) determination of the free energy barrier essential in the predictive use of eq 33. In the following section, we describe a practical method for estimating values for ΔGc to be used in eqs 32 and 33. Two-Dimensional Nucleation Thermodynamics. The free energy barrier ΔGc depends (among other factors) on the shape of a critical nucleus. In this section, we derive closed-form expressions for the minimum free energy shape (Wulff shape) of 2D nuclei found on faces with anisotropic interactions, which can be used in a predictive model to determine the (supersaturation-dependent) steady-state growth shapes of molecular organic crystals via eqs 32 and 33. The free energy associated with the formation of a nucleus from an isothermal−isobaric solution can be expressed in its simplest form as

N i=1

(33)

ΔG(n) ⎛ cnϕe ⎞ 1/2 =⎜ − σn ⎟n kbT ⎝ kbT ⎠

(32)

(36)

where ϕe is an face-specific “effective” edge energy and cn is a shape factor relating n1/2 to the number of molecules located along edges. The shape factor is defined by

While eq 32 is a universal expression for growth by this mechanism, it can only be used in a predictive manner if its constituents can be determined a priori. In eq 32, the nucleation rate J depends on the Wulff shape, whereas the spreading rate depends on f, which depends on a kinetically determined 2D growth shape. The analysis required to predict f has already

cn ≡ 661

P 1/2

A

(37) dx.doi.org/10.1021/cg200855p | Cryst. Growth Des. 2012, 12, 656−669


Article

Figure 3. Free energy barriers for homogeneous 2D nucleation on the (111̅) face of naphthalene grown in ethanol. The PBC network contained in a slice of the face is shown in (a). Possible shapes for 2D nuclei are shown in (b), with the corresponding ΔG(n) shown in (d). Isosurfaces of ΔG(l1,l2,l3) are shown in (c). The Wulff and Wulff-PBC shapes in (b) correspond to the minimum free energy shapes for a given volume unrestricted (12 sides are considered, see the section Does the PBC Approximation Result in an Accurate Prediction of σ2D?) and restricted to PBC directions, respectively, while the kinetic-PBC shape refers to the steady-state growth shape of the nucleus. In (c), the isosurfaces correspond to ΔGc − 0.05, + 0.05, and +1.00kbT for the dark blue, blue, and green surfaces, respectively. In (c) and (d) the solid lines correspond to Wulff (blue) and Wulff-PBC (black) shapes, while the dashed line corresponds to the kinetic-PBC shape. In (d), supersaturations are used for which the Wulff and Wulff-PBC shapes result in ΔGc = 35kbT.

where P is the perimeter of the nucleus with area A. This shape factor will be smallest for disk-shaped nuclei where cn = 2π1/2 ≈ 3.5; for a square nucleus cn = 4. The effective edge energy of a nucleus is given by

ϕe =

Alternatively, a 2D nucleus can be treated as an arbitrary N-sided polygon with the resulting free energy of formation given by N ⎛ edge ⎞ γ ΔG σ = h ∑ ⎜⎜ i ⎟⎟li − H̅ ili kbT kT 2vm i=1 ⎝ b ⎠

∑iN= 1 ni ϕiedge ∑iN= 1 ni

(38)

where H̅ i is the perpendicular distance from the center of the nucleus to edge i (note: the area of an N-sided polygon is given by 1/2∑Ni=1 H̅ ili). The N-dimensional free energy landscape, ΔG(l1,...,lN) generated by eq 41 has a single stationary point, which is a first-order saddle (i.e., there is exactly one negative eigenvalue in the N-dimensional Hessian matrix for ΔG at this point). Such landscapes are shown for six-sided nuclei with 2-fold rotational symmetry (i.e., three pairs of parallel edges) in Figures 3c and 4c. This symmetry allows the reduction from six to three dimensions. The coordinates of the stationary point in the N-dimensional energy surface correspond to the edge lengths and energy barrier for the formation of a critical 2D nucleus. The location of the stationary point can be determined by solving the system of equations given by

where ni is the number of molecules located at the perimeter of the nucleus along edge i and N is the number of edges on the nucleus. The free energy landscape, ΔG(n), along this coordinate (number of molecules) has a single stationary point, which is a maximum (see Figures 3 and 4d). Assuming shapes are constant near this point, there are

nc =

2 2 1 ⎛ cnϕe ⎞ ⎛⎜ 1 ⎞⎟ ⎟ ⎜ 4 ⎝ kbT ⎠ ⎝ σ ⎠

(39)

molecules in a critically sized nucleus. Substituting nc into eq 36, the energy barrier that must be crossed to form a postcritical 2D nucleus is given by 2 ΔGc 1 ⎛ cnϕe ⎞ ⎛⎜ 1 ⎞⎟ = ⎜ ⎟ kbT 4 ⎝ kbT ⎠ ⎝ σ ⎠

(41)

⎞ ⎛∂ ⎜ ΔG(l1, ..., lN )⎟ ⎠l ⎝ ∂li

(40)

Using eq 40 to predict this barrier for 2D nucleation requires the ability to determine the nucleus shape and ϕe a priori.

=0 j≠i

(42)

for i = 1, ..., N. The solution of eq 42 is equivalent to the Wulff construction “discretized” (i.e., restricted) to N-sides, and 662



Article

Figure 4. Free energy barriers for homogeneous 2D nucleation on the (001̅) face of naphthalene grown in ethanol. The PBC network contained in a slice of the face is shown in (a). Possible shapes for 2D nuclei are shown in (b), with the corresponding ΔG(n) shown in (d). Isosurfaces of ΔG(l1,l2,l3) are is shown in (c). In (c), the isosurfaces correspond to ΔGc − 0.05, + 0.05, and +1.00kbT for the dark blue, blue, and green surfaces, respectively. In (d) the blue curves corresponds to the Wulff shape, while the black line corresponds to the Wulff-PBC shape. For this face, the kinetic-PBC shape is indistinguishable from the Wulff-PBC shape. In (d), supersaturations are used for which the Wulff and Wulf-PBC shapes result in ΔGc = 35kbT.

the Frank-Chernov condition, which is given by

is given by

⎛ edge ⎞ v γ H̅ c, i = m ⎜⎜ i ⎟⎟ σ ⎝ kbT ⎠

Gi = const. for i = 1, ..., N Hi

(43)

where Hi is the perpendicular distance from a common point within the crystal to face i.5 For the reasons noted in the section Edge Velocities, it is often impractical to calculate k+ or k− directly, quantities that are necessary for the first-principles prediction of the absolute edge velocities and therefore the prediction of Gi or Hi. However, the Frank-Chernov condition can be modified to

where the subscript c indicates the critical nucleus. A derivation of the unrestricted Wulff shape is given in Appendix A for the case where N is unknown a priori. The critical lengths for each edge on an N-sided nucleus can be determined from eq 43 and the relationship

lc, i =

H̅ c, i + 1 − H̅c , i cos(αi , i + 1)

Ri = const. for i = 1, ..., N xi

sin(αi , i + 1) +

H̅ c, i − 1 − H̅ c, i cos(αi − 1, i) sin(αi − 1, i)

(45)

(46)

where Ri and xi are the relative growth rate and distance from a common point within the crystal to face i, respectively.4 Relative growth rates and distances are defined based on a reference face (denoted by the subscript ref), such that Ri  Gi/Gref and xi  Hi/Href. In doing so, ki+ and ki− are canceled out for each edge (in so much as the assumption of ki+ = k+ for all i is valid; i.e., ki+ = ki for all edges on the crystal) and do not have to be calculated in the shape prediction. Therefore, shapes can be determined using eq 18 to calculate vrel,i (and Rhkl), as opposed to using eq 15 or eq 17 to calculate vi (and Ghkl). However, this approach cannot be used to predict absolute growth rates or the size of a crystal as a function of time. Most of the faces initially present on a crystal “grow out” of existence (becoming edges or vertices4) during the attainment of its steady-state shape. The construction of a crystal’s steadystate shape requires the prediction of the relative growth rates only for the set of faces that are present on this shape. The scaling G ∝ σ2 in the absolute growth rate is the same for

(44)

In order to use eqs 43 and 44 to construct a nucleus, any side determined by eq 44 to have lc,i ≤ 0 is assigned a length of zero, which is a convexity condition. The results of eqs 43 and 44 can be used in eq 41 to determine ΔGc. Practically speaking, the use of this approach requires the a priori determination of the number of edges that should be considered. In the section Does the PBC Approximation Result in an Accurate Prediction of σ2D? we evaluate the errors introduced by restricting the edges initially considered to those coincident with PBC directions. This approach is shown in Figures 3 and 4 for the (111̅) and (001̅) faces of naphthalene grown in ethanol, respectively. Features of these figures are further discussed in later sections. Modeling Multiple Mechanisms. The shape of a growing crystal, with constant perpendicular growth rates for each face, will evolve toward a unique and stable steady state.4,52 The steady-state growth shape of a crystal can be constructed from 663



Article

all F-faces present on the steady-state growth shape when these faces grow exclusively by spirals. In this case, relative growth rates, Ri, are independent of supersaturation. This shape is hereafter referred to as the spiral growth shape. However, G ∝ σ2 is no longer valid for all faces present on the steady-state shape if the supersaturation is high enough that one or more of these faces grows by 2D nucleation. For a face growing by 2D nucleation, the scaling between supersaturation and absolute growth rate is given by

G2D ∝ σ5/6 exp(− F /σ)

approximated by

⎛ σ ⎞1/2 ⎛ Ac, hkl (σ) ⎞1/3 ⎟⎟ ⎜⎜ ⎟⎟ R2D, hkl(σ) ≈ ⎜⎜ ⎝ σ2D, hkl ⎠ ⎝ Ac, hkl (σ2D, hkl) ⎠ ×

where F is defined as

(48)

The scaling of G2D and the definition of F were obtained from eqs 9, 32, 34, 33, and 40 (note: G2D ∝ exp(−ΔGc/3kbT), see also ref 29). The different scalings between G and σ for faces growing by spirals and 2D nucleation results in steady-state shapes that are supersaturation dependent for supersaturations above which at least one face present on the spiral growth shape grows by 2D nucleation. The nature of this dependence agrees with the observation that significant changes to growth shapes often occur over narrow ranges of supersaturations.8−10 Therefore, the ability to account for the different scalings between relative growth rates and supersaturation is necessary in the prediction of supersaturation-dependent growth shapes. In the remainder of this section, we present an approach for doing so that can be rapidly implemented within the framework of a predictive model. As noted in section Growth Mechanisms, at σ2D,hkl the hkl face grows at the same rate by both spiral and 2D nucleation mechanisms; that is, GS,hkl(σ2D,hkl) = G2D,hkl(σ2D,hkl). Using this equality and defining R2D,hkl(σ) as

R2D, hkl(σ) ≡

G2D, hkl(σ) G2D, hkl(σ2D, hkl)

=

G2D, hkl(σ) GS , hkl(σ2D, hkl)

(49)

the absolute growth rates for different faces on the same crystal growing by either a spiral growth or 2D nucleation mechanism can be determined from ⎧ ⎛ ⎞2 ⎪ GS, hkl(σ*)⎜ σ ⎟ for (σ ≤ σ2D, hkl) ⎝ σ* ⎠ ⎪ ⎪ Ghkl(σ) = ⎨ ⎪ R2D, hkl(σ)GS, hkl(σ*) for (σ > σ2D, hkl) ⎪ ⎛ σ2D, hkl ⎞2 ⎟ ⎪ ×⎜ ⎝ σ* ⎠ ⎩

l τS = c c1 v

(52)

(53)

and

τ2D = (Jv 2)1/3 c 2

(54)

for edges with isotropic edge velocities and critical lengths (i.e., vi = v and lc,i = lc for all i), where c1 and c2 are constants relating the coverage and geometries. For low values of σ such that σ ≈ S, the ratio of τS/τ2D is given by

(50)

where σ* is an arbitrarily low supersaturation for which spiral growth is expected for all faces present on the steady-state growth shape. Similarly, the relative growth rate of the hkl face is given by

R hkl(σ) = R2D, hkl(σ)R S, hkl(σ*)

ξhkl(σ2D, hkl)

where Ac is the area covered by a critically sized nucleus. In the application of this model, the requirement of an accurate estimate for σ2D,hkl for each bounding face of the crystal is evident from eqs 49−52. Furthermore, for practical reasons, the reference face should be selected such that σ2D,ref > σ2D for the remaining faces present on the spiral growth shape. Methods for Determining σ2D. Determining σ2D,hkl a priori is equivalent to answering the question: “for what value of σ is GS,hkl = G2D,hkl?” As noted previously, there are several “road-blocks” to estimating absolute growth rates, which make the answer to this question difficult to obtain. Using estimates for edge velocities and face growth rates by a spiral mechanism of v ≈ 5 μm/s and GS ≈ 25 nm/s,5 coupled with the scaling G2D ∼ h(Jv2)1/3, Lovette and Doherty proposed estimated bounds of ΔGc,2D ≲ 35 ± 10kbT for growth by 2D nucleation to occur at a rate comparable to that of spiral growth for molecular organic crystals.29 Furthermore, based on observed transitions between growth by spirals and 2D nucleation, an average of ΔGc,2D = 61 ± 6kbT was determined for the variety of materials in ref 29, none of which were molecular organic crystals. The experimental value is offered to provide perspective for the estimate of ΔGc,2D = 35 ± 10kbT applied herein. In this section, we provide an approximate expression for σ2D (based on several approximations) and determine appropriate ranges of ΔGc,2D for typical molecular organic crystals in order to test the validity of ΔGc,2D ≲ (35 ± 10)kbT. When GS,hkl = G2D,hkl the ratio of τS,hkl/τ2D,hkl equals unity, assuming the same heights for 2D nuclei and spiral edges. Estimates for τS and τ2D for arbitrary polygonal spirals and nuclei are given by

(47)

2 1 ⎛ cnϕe ⎞ F≡ ⎟ ⎜ 12 ⎝ kbT ⎠

ξhkl(σ)

τS = C σ−7/6 exp(− F /σ) τ2D

(55)

where the dimensionless quantity C is defined as

(51)

1/3 ⎛ c1 ⎞⎛ ϕkink ⎞⎛ κ*2D a p2 ⎞ ⎟ ⎟⎟⎜n ̅ C = 2⎜ ⎟⎜⎜ k− ⎟⎠ ⎝ c 2 ⎠⎝ kbT ⎠⎜⎝

where R2D,hkl(σ ≤ σ2D,hkl) = 1. Equation 51 assumes the rate constants k+ and k− are isotropic, and v ̅ and η1 (quantities which are contained in κ2D and cannot be rapidly determined) are independent of supersaturation. Both cn and ϕe, which are contained in F, are face dependent; therefore, R2D(σ) will vary from face to face with σ. Using eqs 32−34, R2D,hkl can be

(56)

with lc, v, and J determined using eqs 20, 9, and 33, respectively, and κ2D *  κ2Dσ−1/2 (such that C is independent of σ). 664



Article

The face-specific values of σ2D can be determined by solving

M σ7/6 = exp(− F /σ)

was used to estimate σ2D,hkl for each face found on the spiral growth shape. Furthermore, the accuracy in our prediction of σ for each face is given by the range of eq 59 for ΔGc = 25 − 45kbT. Does the PBC Approximation Result in an Accurate Estimate of σ2D? The shape factors of critically sized 2D nuclei, cn, were determined using Monte Carlo simulations by ter Horst and Jansens for (001) faces of Kossel crystals with effective edge energies ranging from ϕedge ≈ 0.79 − 5.24kbT.53 The (001) face of a Kossel crystal consists of two orthogonal PBCs of equal strength. Therefore, if nuclei were formed in accordance with the PBCs found in the face, they would be square with cn = 4.0. However, the shape factors determined by ter Horst and Jansens ranged from cn ≈ 3.54 − 3.68 (noting that cn increased with increased ϕedge).53 This discrepancy between the expected and obtained shapes calls into question the accuracy of approximating critical 2D nuclei as polygons bound by edges in PBC directions. In this section, we evaluate the accuracy in predicting σ2D using the PBC approximation to determine nucleus shapes. Nuclei on the (111)̅ and (001)̅ faces of naphthalene grown in ethanol were used to evaluate the accuracy of this approach (in comparison to the “full” Wulff shape). To enable this test, a program was developed with the ability to estimate energies for arbitrary edges on molecular organic crystals in solution, so long as the crystal structure and solubility parameter data for the solvent are known. This program is difficult to implement within a predictive framework, as it requires user input to determine the orientation and number of edges to be evaluated. The program works by first constructing the three-dimensional (3D) network of repeating strong intermolecular interactions found within the crystal using known crystallographic information (i.e., the space group, unit cell parameters, and atomic positions) listed in the Cambridge Structural Database.54 Within this program, the network of intermolecular interactions can be constructed directly using the atomic positions and point charges as inputs to pairwise interaction potentials taken from the generalized amber force field (GAFF)55 or imported from an external program such as HABIT95 or OPiX.56,57 Once the 3D network has been determined, a slice of this volume is taken for the face in question with a specified center point and thickness (typically given by the intermolecular spacing of that face). The solvent accessible surface area (SASA) and intermolecular interaction energies (dispersive, Coulombic, polarization, and total; hereafter denoted as Ed, Ec, Ep, Etot, respectively) of the strong intermolecular interactions found within this slice are then determined, with the SASA obtained using the ShrakeRupley method.58 This slice is then separated into two halves along the direction of the edge in question from a specified center point. These halves are separated such that they are independent of each other; with no overlapping surface contributions or intermolecular interactions. The total intermolecular interaction energies and SASA(s) for each independent half are then determined in the same manner. The SASA is used in lieu of a simple geometric relationship to determine surface areas in an effort to provide an accurate and generalizable approach. The various components of the cohesive energy density of edge,cryst the crystal along the edge i, γ(d,c,p,tot),i are estimated by

(57)

where M ≡ C−1. Equation 57 has several solutions. The solution that yields real results of reasonable magnitude is given by

σ2D = −

6 F 7 W − 6 FM 6/7 −1

(

7

)

(58)

where W−1 is the −1 branch of the Lambert W function (also known as the “product log” function). F and M depend on the specific face being investigated. Equation 58 provides an estimate of the supersaturation where each face transitions from spiral growth to 2D nucleation. While eq 58 provides approximate values of σ2D, the values of κ2D * and k−, which are contained in eq 56, can only be estimated to within a limited degree of certainty without computationally expensive molecular simulations. Using the “typical” values for molecular organic crystals of c1/c2 = 5, ϕkink = 2kbT, n̅ = 10, κ2D * = 6 (1024) cm−2s−1, ap= 6 (10−7) cm, and k− = (105) s−1, results in M ≈ 5 × 10−5. Face-specific values of F can be accurately predicted using eq 48. For most molecular organic crystals, F ≲ 10. As ΔGc = 3(F/σ)kbT, the ability to accurately determine F using the procedure described in section Two-Dimensional Nucleation Thermodynamics enables the prediction of σ2D with the same level of accuracy for which ΔGc,2D = 35 ± 10kbT is valid for the given system. Since M cannot be readily obtained without the use of computationally expensive molecular simulations (due to the inaccuracies in k− and κ2D), the accuracy of this estimate was evaluated over a range of M and F using eq 58 to form the plots shown in Figure 5. The applicability of the estimate ΔGc,2D = 35 ± 10kbT, with reference to the analytical

Figure 5. Plots demonstrating the applicable ranges of F and M where the estimate of ΔGc,2D = 35 ± 10kbT is valid in determining σ2D,hkl. The values of ΔGc,2D/kbT are written adjacent to each curve. The gray regions indicate values of F and M for which the estimate of σ2D,hkl using eq 59 is valid. The behavior near the origin is shown in (b).

solution in eq 58, is indicated by the gray region in the plots. For the case of naphthalene in ethanol, estimated values of M ≈ 5 × 10−5 and F ≈ 0.75 − 2 are well within this region. Therefore, for the case study presented in this article the approximation given by

σ2D, hkl = 3Fhkl /35

edge,cryst

γ(d,c,p,tot), i = 2

(59) 665

energy of bonds broken area created by the two edges

(60)



Article

fitting point charges, an ab initio solution to the electronic structure was determined using the CBS-4M basis set in Gaussian03.64 Using this approach (basis set, point charges, mixing-rules, dielectric, and force field) to calculate the intermolecular interactions, a predicted lattice energy of −17.2 kcal/mol was obtained in HABIT95.56 This predicted lattice energy is within the error of the experimental value reported by NIST of −18 ± 1 kcal/mol.65 In addition, the polarization energies were calculated using the pixel forcefield66 with the default settings used for each of the OPiX programs.57 The strong (based on their contribution to the lattice energy) interactions calculated on this basis and found within slices of the (111̅) and (001̅) faces (shown in Figures 3a and 4a, respectively) are listed in Table 1.

which for each component is given by edge,cryst

γ(d,c,p,tot), i = 2

∑slice E(d,c,p,tot) − ∑halves, i E(d,c,p,tot) ∑ SASA halves, i − SASA slice

(61)

where ∑sliceE(d,c,p,tot) and ∑halves,iE(d,c,p,tot) are the given components of the intermolecular interaction energies summed over slice and over both halves for edge i, respectively. This approach avoids the problem of distinguishing whether a bond is broken on forming the edge by separating the two halves. The approach works for arbitrary edges regardless of their orientation. Once the cohesive energy densities of the solid have been obtained for the direction in question, the edge energy for the solid in solution is approximated using the macroscopic relationship given by

Table 1. Intermolecular Interactions in Naphthalenea

edge

γiedge = γedge,cryst + γsolvent − W ad, i tot, i tot

(62)

solvent where γtot is the cohesive energy density of the solvent on an edge area basis, and Wad,i is the work of adhesion between the two phases. The cohesive energy density of the solvent is determined using solubility parameters. These parameters have been divided by Barton into their dispersive, polar, and hydrogen bonding components.59 The relationship between cohesive energy density and solubility parameters is given by

1/3 2 γsolvent = εVm δz z

edge

+ 2 γedge,cryst γasolvent p, i

δ2h + δ2p

Ed/kbT

Ec/kbT

Ep/kbT

Etot/kbT

0.508 0.597 0.786 0.802

[110] [010] [112] [101]

−6.2 −4.1 −2.2 −0.7

−1.0 −1.7 −0.2 −0.1

−0.3 −0.2 −0.1 0.0

−7.2 −5.9 −2.4 −0.8

On the basis of the values for solubility parameters in ref 59, and the value of ε given by Kaeble60 for ethanol; γdsolvent = 8.28, solvent γasolvent = 15.18, and γtot = 22.52 erg/cm2. Representative slices and the resulting edge energies for the (111̅) and (001̅) faces of naphthalene in ethanol are provided in Supporting Information. These faces contain three PBCs each. The minimum free energy shapes of nuclei on both faces were determined using eqs 43 and 44 and the edge energies listed in Supporting Information. In order to test the accuracy of restricting nuclei to PBC directions for the (111̅) and (001̅) faces of naphthalene, nuclei on these faces were allowed to be bounded by three additional edges that are not parallel to PBCs. The resulting 12-sided nuclei provide more accurate approximations of the Wulff shaped nuclei for each face. These nuclei are denoted as “Wulff” in Figures 3 and 4. Nuclei restricted to edges that are parallel to the PBC directions in each face, and obey eq 43, are denoted as “Wulff-PBC” in Figures 3 and 4. These nuclei are six-sided. Additionally, sixsided nuclei with shapes corresponding to the predicted relative edge velocities in each PBC direction were constructed by setting H̅ c,i/H̅ c,j≠i = vrel,i/vrel,j≠i; the shape of these nuclei are referred to as “kinetic-PBC” shapes. The shape factors, cn, and effective edge energies, ϕe, for these shapes are listed in Table 2. The onset supersaturations for 2D nucleation, σ2D, were determined from the values of cn and ϕe and the assumption that this onset occurs at ΔGc = 35 ± 10kbT (see section Methods for Determining σ2D). Although the Wulff shape resulted in the lowest free energy barrier for the formation of a critically sized nucleus for the (111)̅ and (001)̅ faces of naphthalene in ethanol, the predicted values of σ2D using the “Wulff”, “Wulff-PBC”, and “kineticPBC” shapes for both faces agree to within the uncertainty introduced by approximating ΔGc,2D = 35 ± 10kbT. This is clearly demonstrated by the values of σ2D listed in Table 2.

(64)

where the polar and hydrogen bonding components are grouped as association interactions with59

δa ≡

direction

red cyan mustard green

The colors correspond to lines in Figure 3a and 4a at T = 284 K. To avoid double counting electrostatic contributions, the polarization energies (determined using the pixel force-field)66 were not included in the total interaction energies (determined using the GAFF).55 The column “distance” refers to the distance between centers of mass.

(63)

edge,cryst solvent γd

distance (nm)

a

where ε is a factor accounting for the presence of an interface (typically ε ≈ 0.1)60 and δz is the solubility parameter associated with the component z, with z = d,p,h,a for the dispersive, polar, hydrogen bonding, and associative components, respectively. For this specific set of force-fields (GAFF and OPiX), the work of adhesion along edge i is determined using a geometric mean approach given by W ad, i = 2 γd, i

color

(65) edge,cryst in determining γtot,i , the quantity edge determining Wad,i . Equation 64 is

To avoid double counting edge,cryst γp,i is only used in not ideal for cases where hydrogen bonding interactions are calculated directly within the solid, as more accurate expressions have previously been developed. However, in using these approaches a “consistency” check should be made to determine whether the force-field used to determine the solid-side interactions accurately predicts the lattice energy of the crystal. In our own research, it was determined that for the example system, the approach described in the following paragraph resulted in the most accurate description of the solid-side interactions. For the example of naphthalene in ethanol, the intermolecular interactions were calculated for its monoclinic polymorph (space group P21/a, Cambridge Structural Database refcode NAPHTA10).61 The dispersive, Columbic, and total intermolecular interactions for the crystal were determined using the GAFF force field55 with Kirkwood-Slater mixing rules,62 and restrained electrostatic potential fit point charges63 with a dielectric of 1.2 applied to columbic interactions. Prior to 666



Article

provided, which uses the proposed modifications to calculate interfacial energies that were described in the previous section. These modifications result in a more accurate low supersaturation shape prediction, which is evident by comparing Figure 7 with Figure 8a from ref 6.

Table 2. Predicted Shape Factors, Energies, and Onset Supersaturations for Critical Nuclei for Naphthalene in Ethanola (111̅) Face shape

cn

Wulff Wulff-PBC Kinetic-PBC

3.60 4.07 4.09

ϕe/kbT 0.83 0.90 0.90 (001̅) Face

cnϕe/kbT

σ2D

2.99 3.64 3.70

0.06 ± 0.03 0.10 ± 0.03 0.10 ± 0.03

shape

cn

ϕe/kbT

cnϕe/kbT

σ2D

Wulff Wulff-PBC kinetic-PBC

3.66 3.80 3.80

1.34 1.42 1.42

4.91 5.40 5.40

0.17 ± 0.07 0.21 ± 0.06 0.21 ± 0.06

T = 284 K. The values of σ2D were determined from ΔGc,2D = 35 ± 10kbT. a

Figure 6. Predicted relative growth rates for naphthalene grown in ethanol. The symbols (a), (b), (c), and (d), indicate the {111̅}, {110}, {201̅}, and {001}, families, respectively.

Therefore, on the basis of this analysis, PBC approximations to critical nuclei shapes appear to predict σ2D to within the accuracy required for modeling supersaturation dependent growth shapes. Said differently, restricting nuclei to PBC directions does not introduce an additional error into the model that is beyond the uncertainty in σ 2D , which accompanies the prediction of σ2D from ΔGc,2D = 35 ± 10kbT. Furthermore, when nucleus shapes are restricted to PBC directions there is a broad range of shapes for which ΔGc is within 1kbT of the values determined for the Wulff-PBC shape, including the kinetic shape. This indicates that there are a variety of “acceptable” approximations for the Wulff shape when restricted solely to PBC directions. This is shown by the green surfaces in Figures 3c and 4c. The proximity in free energies between the “Wulff-PBC” and “kinetic-PBC” shapes indicates that there is no additional (rate limiting) energy barrier that must be overcome in the transition from a precritical nucleus to a postcritical nucleus in its kinetic shape. For this system, naphthalene in ethanol, the (001̅) and (111̅) faces used in this analysis contain the strongest and weakest sets, respectively, of intermolecular interactions among faces present on the steady-state spiral growth shape. It can be shown that the values of cnϕe/kbT and σ2D for the “Wulff-PBC” shapes of the remaining faces found on the spiral growth shape (i.e., the {110} and {201}̅ families), are between these limits. Therefore, it is expected that the conclusions from this analysis are valid for all faces present on the steady-state growth shape of the crystal. The accuracy of the PBC approximation to nucleus shapes in determining the energy required to form a critically size nucleus eliminates the need to calculate edge energies in directions which are not coincident with PBCs in the implementation of our model. This enables the rapid prediction of σ2D for each face, without sacrificing an unjustifiable amount of accuracy in comparison to more computationally expensive methods.

Figure 7. Predicted and experimentally obtained steady-state growth shapes for naphthalene in ethanol. The low and mid σ shapes were obtained by Grimbergen et al.10 The scale bar is equal to ≈100 μm for the high σ shape. The flat plate shape in (g) is more clearly observed in the video in WMV format. With the exception of the {001̅} faces in (c), the predictions shown in (c) and (d) should be interpreted as “flat plate” and “spheroid” type shapes rather than shapes bounded by the specific faces shown in the figures. This is because bulk transport becomes the rate limiting process in growth for faces in all but the {001̅} directions for (c) and for all directions in (d).

The following procedure was used to predict supersaturation dependent steady-state growth shapes: 1. The spiral growth shape was constructed using the procedure proposed by Snyder and Doherty,6 which was modified to better account for ϕkink and ϕedge (see eq 64) for a supersaturation of σ = 0.01. 2. For the set of faces present on the spiral growth shape, applying the PBC approximation, 2D nucleus shapes and F were determined using eqs 43, 44, and 48. 3. σ2D,hkl was determined for these faces using eq 59 and a reference face chosen such that σ2D,ref > σ2D for the remaining faces. 4. The values of R2D,hkl(σ) were then determined over the specified range of supersaturations using eq 52. 5. Steady-state growth shapes were constructed at each supersaturation using the values of R2D,hkl determined for each face in eq 51. A bound of Rhkl = 20 was used to estimate the onset of bulk-transport limited growth (red

■

CASE STUDY: NAPHTHALENE GROWN IN ETHANOL Naphthalene grown in ethanol was chosen as a case study to demonstrate the implementation of this model due to the documented effect of supersaturation on growth shapes and its centrosymmetric growth units and PBC network. The nature of this network validates the assumption made in the section Layered Growth of a single type of kink site. Furthermore, a revision of the previously determined6 spiral growth shape is 667



■

APPENDIX A: OBTAINING THE WULFF SHAPE USING CALCULUS OF VARIATIONS The shape corresponding to the lowest free energy at a fixed volume is given by Wulff’s construction, derivations for which are presented by Burton et al.14 and Kern68 among several others. Because of the importance of the result in this article and to clarify an omission/ possible oversight from the earlier work, this construction is derived following a similar method to the aforementioned.14,68 The lowest free energy shape can be obtained by representing the interface of the 2D nucleus as a path, with the result given by eq 73 from ref 68, given by

line in Figure 1). The justification of this estimate is addressed in the Supporting Information accompanying ref 5. The relative growth rates and steady-state growth shapes predicted using this procedure for naphthalene in ethanol are shown as functions of supersaturation in Figures 6 and 7, respectively. The agreement between the predicted steady-state growth shapes and those determined experimentally by Grimbergen et al.10 is shown in Figure 7. Furthermore, although Grimbergen et al.10 did not quantify the supersaturation for which the transition in steady-state shapes occurred (referring instead to the degrees below an unknown solubility at which the solution was held), solubility experiments (performed in our laboratory using the polythermal method)67 indicate that for T = 284 K this transition occurred for σ = 0.09. This is in agreement with the value of σ2D = 0.10 ± 0.03 that was predicted for the {111}̅ family of faces (compare the predictions in 7a and b to the experimental results in Figure 7e,f). Using the apparatus described by Lovette et al.,52 a crystal with a flat-plate shape was grown at a high supersaturation; this crystal is shown in Figure 7g. This shape is similar to the growth shape predicted for σ ≈ 0.14−20, shown in Figure 7c. At this level of supersaturation, bulk transport limited growth was predicted for all faces that are not included in the {001̅} family, based on the upper bound stated in step 5. This prediction is supported by the transition from stable to unstable growth (demonstrated by the loss of convexity in the crystal shape) at this high supersaturation, which is shown in the timelapse video. Although a spherulitic shape with all faces growing at a bulk transport limited rate was predicted for σ > 0.27, this supersaturation is beyond the MZW determined for this system, and therefore, is difficult to obtain in a traditional cooling crystallizer.

■

Article

AP = h

∫0

2π

⎛ k Tσ ⎞ VR = h⎜ b ⎟ ⎝ 2vm ⎠

γ(θ)(H̅ (θ) + Ḧ (θ)) dθ

∫0

2π

H̅ (θ)(H̅ (θ) + Ḧ (θ)) dθ (66)

where H̅ (θ) and Ḧ (θ) are the distances and second derivatives of distance with respect to θ of the 2D nucleus from its center, respectively, and γ(θ) is the edge energy at an inclination θ, as measured from a reference angle. Calculus of variations can be used to determine the lowest free energy with the constraint of a fixed volume. Defining

P(γ(θ), H̅ (θ), Ḧ (θ)) ≡ h γ(θ)(H̅ (θ) + Ḧ (θ)) Q (H̅ (θ), Ḧ (θ)) ≡ h(kbT σ/2vm)H̅ (θ)(H̅ (θ) + Ḧ (θ))

(67)

the resulting Euler−Lagrange equation is given by

⎛ ∂Q ⎞ ∂P ⎟=0 P − Ḧ + λ⎜Q − Ḧ ⎝ ∂Ḧ ∂Ḧ ⎠

(68)

where λ is the Lagrange multiplier. Solving eq 68 for H̅ (θ) yields

SUMMARY

H̅ (θ) =

In this article, we have provided a means for incorporating discrete supersaturation effects into a rapidly implementable predictive model for steady-state growth shapes. In doing so, a generalized model for 2D nucleation was proposed that can be integrated with previous predictive models for spiral growth shapes.6 Furthermore, the assumptions of critically sized nuclei being bounded solely by PBC directions and the prediction of σ2D on the basis of ΔGc,2D = 35 ± 10 kbT are carefully justified. These assumptions enable the rapid implementation of this model. The efficacy of this model for predicting the influence of supersaturation on steady-state growth shapes was demonstrated in a case study of naphthalene grown in ethanol. In applying this model to systems with non-centrosymmetric PBCs,18 further modifications are necessary to account for the possibility of different slice energies between successive layers and edges aligned in the same direction with different energies. Furthermore, in light of the recognition that the Voronkov velocity profile (given by eq 23) is a poor estimate for edges with ϕkink ≲ 5kbT, modifying the closed form expression for rotai tion times applied in this model to determine RS that assumes consistency with this profile is a potential area for further improvement.

2 vm ⎛ γ(θ) ⎞ ⎟ ⎜ λ σ ⎝ kbT ⎠

(69)

Lagrange multiplier of the critical nucleus can be determined through the substitution of H̅ (θ) into eq 66, yielding

1 1 − APc + APc = APc (70) λ 2 which results in λ = 2 (note: ΔGc = 1/2APc is well-known14). Therefore, for a critically sized nucleus in its lowest energy shape

H̅ (θ) =

vm ⎛ γ(θ) ⎞ ⎟ ⎜ σ ⎝ kbT ⎠

(71)

Equation 71 is equivalent to the Wulff construction upon the introduction of a further constraint ensuring convexity, given by

■

γ(θ) + γ̈(θ) ≥ 0

(72)

ASSOCIATED CONTENT

* Supporting Information S

Diagrams showing representative slices of the (111)̅ and (001)̅ faces of naphthalene in ethanol, and a table containing the corresponding edge energies. This information is available free of charge via the Internet at http://pubs.acs.org/. 668



Article

W Web-Enhanced Feature *

(30) Teng, H. H.; Dove, P. M.; Orme, C. A.; De Yoreo, J. J. Science 1998, 282, 724−727. (31) Chernov, A. A.; De Yoreo, J.; Rashkovich, L. N.; Vekilov, P. G. MRS Bull. 2004, 29, 927−934. (32) Voronkov, V. V. Sov. Phys. Cryst. 1973, 18, 19−23. (33) Liu, X. Y.; Bennema, P. Phys. Rev. B 1996, 53, 2314. (34) Mighell, A. D.; Rodgers, J. R. Acta Crystallogr., Sect. A 1980, 36, 321−326. (35) Dunitz, J. D.; Gavezzotti, A. Chem. Soc. Rev. 2009, 38, 2622− 2633. (36) Nielsen, A. E. Kinetics of Precipitation; Pergamon Press Limited: Oxford, 1964. (37) Lewis, B. J. Cryst. Growth 1974, 21, 29−39. (38) Cuppen, H.; Meekes, H.; van Enckevort, W.; Vlieg, E. J. Cryst. Growth 2006, 286, 188−196. (39) Gilmer, G. H.; Bennema, P. J. Appl. Phys. 1972, 43, 1347−1360. (40) Gilmer, G. H. J. Cryst. Growth 1976, 36, 15−28. (41) Cuppen, H. M.; van Veenendaal, E.; van Suchtelen, J.; van Enckevort, W. J. P.; Vlieg, E. J. Cryst. Growth 2000, 219, 165−175. (42) Rak, M.; Izdebski, M.; Brozi, A. Comput. Phys. Commun. 2001, 138, 250−263. (43) Boerrigter, S.; Cuppen, H.; Ristic, R.; Sherwood, J.; Bennema, P.; Meekes, H. Cryst. Growth Des. 2002, 2, 357−361. (44) Boerrigter, S.; Josten, G.; vandeStreek, J.; Hollander, F.; Los, J.; Cuppen, H.; Bennema, P.; Meekes, H. J. Phys. Chem. A 2004, 108, 5894−5902. (45) Piana, S.; Reyhani, M.; Gale, J. D. Nature 2005, 438, 70−73. (46) Deij, M. A.; Cuppen, H. M.; Meekes, H.; Vlieg, E. Cryst. Growth Des. 2007, 7, 1936−1942. (47) Falo, F.; Bishop, A. R.; Lomdahl, P. S.; Horovitz, B. Phys. Rev. B 1991, 43, 8081. (48) Piana, S.; Gale, J. D. J. Am. Chem. Soc. 2005, 127, 1975−1982. (49) Piana, S.; Jones, F.; Gale, J. D. J. Am. Chem. Soc. 2006, 128, 13568−13574. (50) Szurgot, M.; Prywer, J. Cryst. Res. Technol. 1991, 26, 147−153. (51) Kashchiev, D. Nucleation: Basic Theory with Applications; Butterworth-Heinemann: Oxford, 2000. (52) Lovette, M. A.; Muratore, M.; Doherty, M. F. AIChE J. 2011, No. DOI: 10.1002/aic.12707. (53) ter Horst, J.; Jansens, P. Surf. Sci. 2005, 574, 77−88. (54) Allen, F. H. Acta Crystallogr., Sect. B. 2002, 58, 380−388. (55) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (56) Clydesdale, G.; Roberts, K. J.; Docherty, R. J. Cryst. Growth 1996, 166, 78−83. (57) Gavezzotti, A. OPiX, A Computer Program Package for the Calculation of Intermolecular Interactions and Crystal Energies; University of Milano, 2003. (58) Shrake, A.; Rupley, J. A. J. Mol. Biol. 1973, 79, 351−364. (59) Barton, A. F. M. Chem. Rev. 1975, 75, 731−753. (60) Kaelble, D. H. Physical Chemistry of Adhesion; WileyInterscience: New York, 1971. (61) Brock, C. P.; Dunitz, J. D. Acta Crystallogr., Sect. B 1982, 38, 2218−2228. (62) Leach, A. R. Molecular Modelling Principles and Applications, 2nd ed.; Pearson Education Limited: Harlow, England, 2001. (63) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (64) Frisch, M. J. et al. Gaussian 03, Revision D.02; Gaussian, Inc.: Wallingford, CT, 2004. (65) Linstrom, P., Mallard, W., Eds. NIST Chemistry WebBook, NIST Standard Reference Database Number 69; National Institute of Standards and Technology, Gaithersburg MD, http://webbook.nist. gov, November 14, 2009. (66) Gavezzotti, A. J. Phys. Chem. B 2002, 106, 4145−4154. (67) Fowler, D. Unpublished. (68) Kern, R. In Morphology of Crystals: Part A; Sunagawa, I., Ed.; Terra Scientific Publishing Company: Tokyo, 1987; Chapter 2, pp 77−206.

A video in WMV format of naphthalene growth in ethanol is available in the HTML version of the paper.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

■

ACKNOWLEDGMENTS The authors wish to acknowledge the financial support of Eli Lilly. The authors thank David Fowler for the micrograph shown in Figure 7g.

■

REFERENCES

(1) Winn, D.; Doherty, M. F. AIChE J. 1998, 44, 2501−2514. (2) Winn, D.; Doherty, M. F. AIChE J. 2000, 46, 1348−1367. (3) Rohl, A. L. Curr. Opin. Solid State Mater. Sci. 2003, 7, 21−26. (4) Zhang, Y.; Sizemore, J. P.; Doherty, M. F. AIChE J. 2006, 52, 1906−1915. (5) Lovette, M. A.; Browning, A. R.; Griffin, D. W.; Sizemore, J. P.; Snyder, R. C.; Doherty, M. F. Ind. Eng. Chem. Res. 2008, 47, 9812− 9833. (6) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. A 2009, 465, 1145− 1171. (7) Hartman, P.; Perdok, W. G. Acta Crystallogr. 1955, 8, 49−52. (8) Durbin, S.; Feher, G. J. Cryst. Growth 1991, 110, 41−51. (9) Shekunov, B.; Grant, D. J. Phys. Chem. B 1997, 101, 3973−3979. (10) Grimbergen, R.; Reedijk, M.; Meekes, H.; Bennema, P. J. Phys. Chem. B 1998, 102, 2646−2653. (11) Ristic, R.; Finnie, S.; Sheen, D.; Sherwood, J. J. Phys. Chem. B 2001, 105, 9057−9066. (12) Land, T. A.; De Yoreo, J. J. J. Cryst. Growth 2000, 208, 623−637. (13) Sleutel, M.; Maes, D.; Wyns, L.; Willaert, R. Cryst. Growth Des. 2008, 8, 4409−4414. (14) Burton, W. K.; Cabrera, N.; Frank, F. C. Phil. Trans. R. Soc. A 1951, 243, 299−358. (15) Voronkov, V. Sov. Phys. - Crystallogr. 1970, 15, 8−13. (16) Markov, I. V. Crystal Growth for Beginners, Fundamentals of Nucleation, Crystal Growth and Epitaxy; World Scientific: Singapore, 2003. (17) Vekilov, P. G. Cryst. Growth Des. 2007, 7, 2796−2810. (18) Kuvadia, Z. B.; Doherty, M. F. Cryst. Growth Des. 2011, 11, 2780−2802. (19) Lovette, M. A.; Doherty, M. F. Phys. Rev. E 2011, in press. (20) Sizemore, J. P.; Doherty, M. F. J. Cryst. Growth 2010, 312, 785− 792. (21) Chernov, A.; Komatsu, H. In Science and Technology of Crystal Growth; van der Eerden, J., Bruinsma, O., Eds.; Kluwer Academic Publishers: Boston, 1995. (22) Ohara, M.; Reid, R. C. Modeling Crystal Growth Rates from Solution; Prentice-Hall, Inc.: Upper Saddle River, NJ, 1973. (23) Chen, K.; Vekilov, P. G. Phys. Rev. E 2002, 66, 021606. (24) Petsev, D. N.; Chen, K.; Gliko, O.; Vekilov, P. G. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 792−796. (25) Vekilov, P. In Nanoscale Structure and Assembly at Solid-Fluid Interfaces; Liu, X. Y., De Yoreo, J. J., Eds.; Kluwer Academic Publishers: Boston, 2004; Vol. II: Assembly in Hybrid and Biological Systems; pp 145−200. (26) Vekilov, P. In Protein Nanotechnology: Protocols, Instrumentation and Applications; Vo-dinh, T., Ed.; Humana Press: New York, 2005; pp 15−52. (27) Hänggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 1990, 62, 251. (28) Paloczi, G. T.; Smith, B. L.; Hansma, P. K.; Walters, D. A.; Wendman, M. A. Appl. Phys. Lett. 1998, 73, 1658−1660. (29) Lovette, M. A.; Doherty, M. F. J. Cryst. Growth 2011, 327, 117− 126. 669


Predictive Modeling of Supersaturation-Dependent Crystal Shapes

Recommend Documents