Crystal

A Structural-Kinetic Approach to Model Face-Specific Solution/ Crystal Surface Energy Associated with the Crystallization of Acetyl Salicylic Acid from Supersaturated Aqueous/Ethanol Solution Robert B. Hammond, Klimentina Pencheva, and Kevin J. Roberts*

CRYSTAL GROWTH & DESIGN 2006 VOL. 6, NO. 6 1324-1334

Institute of Particle Science and Engineering, School of Process, EnVironmental and Materials Engineering, UniVersity of Leeds, Leeds, LS2 9JT, UK ReceiVed October 20, 2005; ReVised Manuscript ReceiVed February 22, 2006

ABSTRACT: Classical homogeneous nucleation theory has been integrated with molecular modeling techniques to model the effects of cluster shape-anisotropy associated with nucleation from solution. In this approach, the geometric shape of a crystal nucleus is modeled assuming it equates to the predicted growth morphology of the resultant macroscopic crystal with the later simulated via an attachment energy model using empirical intermolecular force calculations adopting the atom-atom approximation. A new coupled model integrating nucleation theory, morphological simulation and solvent binding calculations, has been developed and combined with experimental nucleation data for the case example of acetyl salicylic acid (aspirin) crystallizing from an aqueous/ ethanol solution. Nucleation parameters, such as critical nucleus size and specific surface energy, have been calculated and compared for the cases of both isotropic (spherical) and anisotropic (polyhedral) nucleus models. A comparison between nucleation data modeled on the basis of both shaped polyhedral and spherical clusters reveals a larger critical nucleus size for the anisotropic (Dequivalent ) 16.54 Å at 40 °C) compared to the isotropic (D ) 13.07 Å at 40 °C) shape model. Theoretical molecular modeling calculations of the specific surface energy anisotropy in solution show the specific surface energy for the habit planes of aspirin to vary from 3.65 mJ/m2 for the dominant {100} facet to 113 mJ/m2 for the minor {1h11} facet. A comparison between the dominant crystal habit planes, notably, the hydrophilic {100} and hydrophobic {002}, reveals them to be more differentiated in terms of their specific surface energy in solution (3.65 and 10.90 mJ/m2) than in a vacuum (829 and 904 mJ/m2), respectively, in good agreement with their known surface chemistry. The total, specific surface energy (41.9 mJ/m2) in solution calculated via molecular modeling for the polyhedral-shaped clusters was found to be somewhat larger, but still in pleasing general agreement with that calculated from experimental nucleation data as determined using induction time measurements assuming a spherical nucleus (3.99 mJ/m2). The potential for the further development of this overall modeling approach is reviewed. 1. Introduction An important step in the production of a new solid phase from solution via crystallization is the formation of metastable molecular clusters of a critical size, i.e., nucleation, as a precursor to crystal growth. Controlling the relative importance of the nucleation and growth processes in crystallization unit operations is often a key issue in the scale-up of material systems from laboratory through to manufacturing scale sizes. Hence, there is currently wide-ranging interest in understanding the molecular-scale details of nucleation processes. There are two types of nucleation processes associated with formation of a crystalline phase from solution or melt: first, homogeneous nucleation, production of molecular clusters of increasing size in the bulk of the mother phase that have structures similar to the structure of the emerging crystalline phase and second, heterogeneous nucleation, production of molecular clusters on the surfaces of impurity species so reducing the free energy of activation associated with formation of a critical cluster. However, experimental discrimination between homogeneous and heterogeneous nucleation mechanisms remains challenging. The description of homogeneous nucleation using fluctuation theory is long standing,1 but, in particular, potential applications dominated by the occurrence of homogeneous nucleation, as opposed to heterogeneous nucleation, is small.2,3 Despite this, experimental studies applying classical homogeneous nucleation theory4-6 are still useful as they provide comparability between different crystallizing environments, e.g., different chemical systems, process scale-sizes, and the influence of impurities. * To whom correspondence should be addressed. Tel: +44(0)113 3432408. Fax: +44(0)113 3432405. E-mail: [email protected].

The theory underpinning classical homogeneous nucleation has been a subject of continuous study and debate during the past 20 years7-13 reflecting the inherent difficulty in validating the theory due to problems associated with measuring the sizes of nuclei. As a result, the paucity of reliable nucleation data, and hence measurements of the specific surface energy at the crystal nuclei/solution boundary, is long standing. Therefore, our ability to predict the number of molecules within nucleation clusters based on existing theory remains quite limited. A further issue relates to the lack of a theoretical approach that takes into account the expected structure of nucleation clusters with respect to both the anisotropy of the intermolecular forces and the potential for nucleation clusters to provide a molecular template for the formation of the crystalline phase post-nucleation. In particular, the inter-relationship between molecular-scale nuclei structure and variation in the crystal phase (polymorphism) remains unclear at this time and requires further detailed study at a fundamental level. Set within the above context, this paper seeks to develop and apply a more molecularly based model of the nucleation process through the integration of classical nucleation theory with molecular modeling procedures associated with morphological prediction. Hence, the starting point for this study draws on an assumption that instead of adopting an isotropic (spherically shaped) cluster structure, the shape of the nucleus could be better modeled using a material’s final crystal growth morphology based on the crystallography of the subsequent solid-state phase following growth. This combined nucleation/structure model has been developed and confronted with measured nucleation kinetic data associated with the crystallization of acetyl salicylic acid (aspirin) from solution in an aqueous-ethanol solution.

10.1021/cg0505618 CCC: $33.50 © 2006 American Chemical Society Published on Web 04/25/2006

Modeling of Face-Specific Solution/Crystal Surface Energy

Crystal Growth & Design, Vol. 6, No. 6, 2006 1325

2. Background Theory 2.1. Nucleation and Nucleation Rate. For solution growth, the formation of a crystalline phase involves aggregation of the solute molecules to form clusters driven by solution supersaturation. As soon as such clusters grow to a critical size, bulk intermolecular forces among the molecules within the cluster begin to prevail over the cluster’s surface interactions with molecules in the surrounding crystallization solution creating a viable form for subsequent crystal growth.14 The number of molecules expected to be in such a stable nucleus depends on a number of process and system-specific parameters and has been estimated to range from about 10 to several thousand molecules.15 When a cluster of volume V and surface area S is formed, the associated change in the cluster free energy (∆G) is given by

V ∆G ) GV + GS ) - kBT ln σ + Sγ ν

(1)

where ν is molecular volume, kB is the Boltzman constant, T is the absolute temperature, σ is the supersaturation defined as a ratio between the actual concentration and saturation concentration, and γ is the solution specific surface energy or interfacial tension. (Note: As this work focuses on crystallization from solution and hereafter the term solution (specific) surface energy is used.) The first term in eq 1 expresses the difference in the chemical potential between the crystalline phase and the surrounding mother liquor and the second term describes the free energy of formation of the new interface. So the free energy of the system decreases in proportion to the volume of new phase created but increases in proportion to the area of the interface formed. Hence, for a given supersaturation, the formation of a viable nucleus of the new phase not only requires overcoming a certain free energy barrier but also, in doing so, realizing a particular molecular cluster configuration from among the enormous number that are plausible. In the classical approach, the nucleus cluster (see, for example, ref 16) is considered to be spherical in shape and so, after substituting the volume and the surface area for a sphere, differentiating eq 1 with respect to the cluster size (D) and equating the outcome to zero:

d∆G )0 dD

(2)

an equation for the critical free energy change (∆Gcritical), associated with the formation of a cluster with a spherical shape and size Dcritical can be derived (Figure 1). (Note: To provide clear linkage between nucleation cluster size and experimental assessment of particle size cluster diameter rather than radius have been standardized in terms of terminology.)

πγDcritical2 ∆Gcritical ) 3 Dcritical )

4γν kT ln σ

Figure 1. Schematic diagram showing the free energy balance between the surface area and volume terms with the maximum corresponding to a critical nucleus size (Dcritical).

(3) (4)

The rate of nucleation can be defined as the number of nuclei formed per unit time per unit volume and using the Arrhenius rate equation together with the Gibbs-Thomson

relationship15 an expression for the nucleation rate (J) can be derived

(

)

16πγ3ν2 J ) A0 exp - 3 3 3kB T (ln σ)2

(5)

where A0 is the preexponential factor relating to the frequency of inter-molecular collision. Hence, the nucleation rate is governed by three main variables: temperature, degree of supersaturation, and the solution specific surface energy. 2.2. Crystal Size and Morphology. A limitation of the use of classical nucleation theory for calculating the specific surface energy reflects the fact that it relies on the assumption that the critical nucleus can be modeled using a simple geometric model, i.e., spherical shape. While this restriction can be mitigated to a degree through the use of empirical shape factors, these do not rigorously take into account the structural chemistry associated with the formation of nucleation clusters, relying instead upon a simple model based on a geometrical approximation. Mindful of the shape anisotropy for crystallographic materials in general and in particular those of lower crystal symmetry, it is important to rationalize the concept of the crystal “size” within the framework of a polyhedral shaped nucleus. For a smooth spherical nucleus, it is obvious that there is only one characteristic dimension, i.e., the nucleus diameter. However, for an irregular object, such as the polyhedron, this definition becomes less clear and the approach usually adopted is to use either an equivalent volume (Dv) or surface diameter (Ds)17 given by

Dv )

x6Vπ

(6)

Ds )

xπS

(7)

3

2

where V is the particle volume and S is its surface area. Previous studies have examined nonspherical nuclei clusters, e.g., Chenthenmari18 assigned a hyperbolical shape to the nuclei to favor antiparallel alignment of the dipoles to weaken dipoledipole interactions. Van Gelder19 used a cylindrical nucleus model to describe the nucleation of long-chain alkyl surfactant salts. More usually though the nucleus is neither a sphere nor some other regular shape but likely to be a cluster of molecules

1326 Crystal Growth & Design, Vol. 6, No. 6, 2006

Hammond et al.

having an irregular polyhedral defined shape. The latter can be modeled based upon the assumption that crystallizing nuclei would be likely to adopt a particle shape consistent with the nature of its intermolecular bonding motif. In this respect, a facetted polyhedral model of a nucleus, mimicking the shape of the developed crystal, is attractive in that it provides, potentially, a molecularly-based method for examining crystal nucleation in a self-consistent manner. The basic procedures for the prediction of crystal morphology of molecular crystals based on their crystallographic structures are now well established.20 In this, the relative growth rate (Rgrowth) of the crystal can be taken to be proportional to the surface attachment energy Eatt using the Hartman-Bennema21 approximation, i.e.,

calculate the surface energies using a bond-breaking model.32 Thermodynamic considerations have led to a proportionality relationship between the cluster size and the surface energy in a manner which takes into consideration the significant contribution made by the crystal corners and edges particularly at small particle sizes.32 However, as its variation decreases quite rapidly with the increase of the crystal size,32,33 surface energy is assumed size-independent in the classical nucleation theory.8 Numerical values of surface energies can be estimated using attachment energy calculations21 thus,

Rgrowth ∝ Eatt ) Elatt - Esl

where Z is the number of molecules in the unit cell of volume Vcell, Eatt is the attachment energy (for definition of this parameter see section 2.3), dhkl is the d spacing, NA is Avogadro’s number. It should be noted that the relationship in eq 9 is only strictly speaking true for crystallographic forms having large d spacings reflecting the situation described here as morphologically important faces usually have large d spacings. Another common method for specific surface energy determination is to use experimentally determined induction time data5,34 based on the assumption that the induction time can be taken as a measure of the inverse of the nucleation rate.

(8)

where Eatt is defined22 as the energy released on addition of the growth slice (dhkl) to the growing crystal surface, Elatt is the crystal lattice energy given as the sum of all intermolecular interactions within the solid state, and Esl, the slice energy, is the fraction of these interactions that lie within the most elementary lattice plane spacing (dhkl). The lattice plane spacing can be derived using the classical Bravais-Friedel-DonnayHarker (BFDH) rules (see, e.g., refs 23 and 24) in which the effect of translational lattice symmetry is taken into account to identify the fundamental surface layer growth systems. From this, the predicted polyhedral form of the crystal can be computed via a polar plot of the Eatt(hkl) made as a function of crystal orientation using the gnomonic projection25 with the overall methodology providing good general agreement with observed crystal morphologies. A number of studies applying this approach to a range of morphological predictions for organic compounds have been reported.20,24,26-28 2.3. Surface Energy Determination. The concept of surface energy is best explained by drawing a distinction with the volume energy associated with molecules that reside within the bulk of a crystal lattice and that are subjected to forces exerted by surrounding molecules in 3D. In contrast, for a molecule at a crystal surface in equilibrium with a vacuum the forces are rather anisotropic in nature, and hence, for the surface molecules to be at equilibrium, a different intermolecular arrangement results at the surface when compared to that in the bulk. This leads to these molecules possessing additional (surface) energy in comparison to the corresponding molecules inside the bulk solid.29 This total surface energy will be proportional to the surface area, whereas the specific surface energy is related to unit surface area. The latter will depend on the packing arrangement of the molecules within the surface for the particular face exposed at the surface. Since the surface molecular arrangement differs as a function of crystal orientation (hkl), the specific surface energy of different crystal habit-faces will, in turn, have different values. When a crystalline particle is surrounded by wetting mother liquor, the surface molecules will, in this case, differ from those in the vacuum case, as described above, because they are also attracted by liquid-phase molecules, hence reducing the energy of the intermolecular interactions between surface active molecules. This has the effect of decreasing the surface energy, e.g., the specific surface energies for the major faces of CaCO3 have been calculated30 to range from 1089 to 1787 mJ/m2, while experimental values for solution crystallization are considerably lower, e.g., 100 mJ/m2. Some predictions of the specific surface energy have been made by applying molecular dynamic calculations,31 while Walton has reported a simple method to

γhkl )

ZEattdhkl 2VcellNA

(9)

3. Experimental and Computational Methods The overall methodology developed and applied in this work introduces two major components: (a) Experimental determination of the critical size and specific surface energy of an isotropic (spherical) nucleus. (b) The anisotropic (polyhedral) nucleation model containing three sub-components: (i) simulation of the crystal morphology; (ii) calculation of the surface area and volume for a polyhedral nucleus; (iii) calculation of surface anisotropy for the specific surface energy of facetted crystal nuclei. This overall modeling methodology is summarized in Figure 2 and a detailed description of all the stages is provided in this section. 3.1. Experimental Nucleation Studies Associated with Calculation of Specific Surface Energy. The material studied was acetyl salicylic acid (aspirin) with a molecular weight of 180.2 and purity of 99.0% recrystallized from ethanol with a molecular weight of 46.07 and purity of 99.8% both purchased from Sigma Aldrich Company. The solubility of aspirin in the mixed water/ethanol solvent was taken from literature data.35 Aspirin was recrystallized from 38% ethanol-water solutions under isothermal conditions, and induction time measurements were carried out at different temperatures. The experiments were carried out using a HEL (www.helgroup.co.uk) 100 mL AUTOMATE jacketed crystallizer equipped with temperature and turbidity probes. After the material was completely dissolved in the solvent, the system was crash cooled with a cooling rate of 1.3°/min to a predetermined temperature and held under isothermal conditions until the turbidometric detection of the appearance of the first nuclei was observed at which point the elapsed time period (τ) was noted. The induction time was taken to be proportional to the inverse of the nucleation rate using the homogeneous rate equation (eq 5), thus:



Figure 2. Schematic diagram highlighting the overall methodology associated with the polyhedral nucleation model and POLYPACK program together with a summary of the main interaction pathways for the different sub components used.

ln τ ) ln B0 +

16πγ3ν2NA3 3R3T3(ln σ)2

(10)

Using this and plotting the logarithm of τ versus 1/T3(ln σ)2 provides a linear plot from which the specific surface energy can be estimated from a slope of the line and the preexponential factor from the intercept. The size of the critical nucleus was determined as a function of temperature and solution supersaturation through eq 4. The shape of the crystals obtained in these experiments was assessed using an inverted Olympus IMT-2 microscope equipped with a video camera and image analysis system. 3.2. Morphology Prediction and Calculation of Surface Energy. The crystal morphology was predicted using the attachment energy method20,23 in which the intermolecular interactions contributing to the crystal lattice energy (Elatt) were partitioned into surface slice (Esl) and attachment (Eatt) energy. This was carried out using the HABIT98 program,36 taking as input the crystallographic structure together with selection of low index habit planes (hkl) as derived from the BFDH rules. The intermolecular forces were calculated using the atomatom method together with an appropriate potential function, in this case derived by Momany et al.37 force field, which was validated via a comparison of the resultant calculated lattice energy with available sublimation enthalpy data thus:

Elatt ) ∆Hsub - 2RT

(11)

In these calculations, Ewald summation was not needed reflecting the small electrostatic contribution to the lattice sum.38 Attachment energies for the (hkl) surfaces selected via the BFDH rule were calculated based on minimizing the slice energy with respect to the slice boundaries to ensure that the most stable (slowest growing) slice was modeled. The surface relaxations were not explicitly taken into account for the calculation of attachment and subsequently surface energy as such effects are not as significant for organic39 as compared with inorganic40 compounds. Identification of the structural unit that adds to growing crystal surfaces requires careful consideration of both the molecular structure of the solute material and the solvent environment. A number of authors41-43 have shown, however, that solutions of polar solvents are more likely to contain solute molecules as

monomers rather then dimers as is the case for the system studied here. Hence in these calculations, it was assumed that the aspirin molecules in solution are dispersed as monomeric species reflecting the polar solvent environment and dock accordingly. 3.3. Polyhedral Volume and Surface Area as a Function of Size. A dedicated FORTRAN program POLYPACK (full details are given in Supporting Information) was used to calculate the polyhedral shape based on inputted attachment energies. In this, matrix algebra calculations,44 based on the methodology outlined by Dowty,45 were used to calculate the Cartesian coordinates of the corners of the smallest possible irregular polyhedron enclosed by specified (hkl) planes. The volume and the surface area of this polyhedron were modeled as a function of crystal size providing parameters needed for the nucleation prediction. The size-specific function was achieved via multiplication of the polyhedral coordinates using an expansion coefficient based on a root morphological simulation at unit cell size scale. Hence, the smallest possible polyhedron enclosed by the dominant morphological faces (hkl), in which the perpendicular center to face distances are taken to be proportional to the attachment energies, has been expanded to represent different polyhedral sizes. In this, the expansion coefficient simply scales all linear dimensions for the predicted polyhedral form, thus producing a polyhedron of any specific dimensions. Using this, polyhedral volume and surface area as a function of this expansion coefficient can easily be calculated. Figure 3 shows an example of its use in the calculation of polyhedral surface and volume equivalent diameters using eqs 6 and 7 for the case of aspirin. This enabled, using an appropriate coefficient, a particular size of the polyhedron to be obtained. In the procedure adopted, the irregular polyhedron was modeled with respect to two reference spheres, one calculated to have the same volume and the other calculated to have the same surface area. The volume and surface energies for different sizes of polyhedral nuclei were determined using of eq 1 together with the modeled polyhedral volume and surface areas. Using data obtained from these calculations, the Gibbs free energy as a function of the crystal size, and hence the critical cluster size, of the polyhedral nucleus was determined. 3.4. Anisotropic Specific Surface Energy as a Function of Habit Face (hkl). The anisotropy of the solution specific surface energy was directly modeled via a partition of the total


Hammond et al.

Figure 3. Calculated polyhedral expansion coefficient as a function of equivalent diameter using volume and surface area equivalent models for aspirin crystal.

surface energy between the individual crystal habit faces. The methodology for this, shown schematically in Figure 4, involves calculation of the surface areas of the individual crystal habit (hkl) faces (Si) and total surface area (Stotal), which can be directly related to the specific surface energy (γ), considered as an appropriate average11 of the individual specific surface energies (γi) for the crystal faces as presented in the predicted crystal morphology thus N

Stotalγ )

∑1 MiSiγi

(12)

where M is the multiplicity of the crystal face and N is the number of forms displayed in the crystal morphology. One has to bear in mind that an approximation of equating the surface energy (the potential energy the surface stores) and surface free energy (energy released for formation of new surface) has been made. The surface energy given by

Es ) Gs + TSs

(13)

is the sum of the surface free energy and the surface entropy multiplied by the absolute temperature. An assumption of this work is that in terms of surface binding calculations the entropy term can be neglected. This assumption is probably acceptable providing the growth temperatures are comparatively low and the potential for a wide range of molecular conformation states for the surface adsorbed species is not significant. The specific surface energy for the individual faces can be calculated easily based upon a vacuum level molecular modeling simulation using eq 9. However, in the real case, the specific solution surface energy will be decreased due to the wetting of the crystal by the solution, and hence for solution growth the interactions with the solvent need to be explicitly taken into account. To assess surface solvent wetting, molecular modeling using Cerius2 software46 and a surface docking approach47,48 were used to estimate the solvent binding energy for the individual habit (hkl) crystal faces. In this, a molecular model of the surfaces was built up from the bulk crystallographic structure as three layers of molecules of which the top one was subjected to surface relaxation. In contrast to attachment energy calculations

(section 3.2), surface relaxation effects can be expected to be quite important in assessing the energy of an isolated solvent or host molecule binding to the crystal surface.49 For all plausible binding configurations to be taken into account and to find the global minimum of the system, a number of starting configurations of a docked molecule on the crystal surface were generated with the aid of the SYSTSEARCH program50 using a grid-based search method.51 In this procedure, the solvent/host molecule was placed into position on the cleaved crystal surface, having a minimized top layer (see Figure 5), and translated in two dimensions, with a step size equal to 5 Å for the host aspirin molecule, and 2 Å for the two solvent molecules and rotated with an angular step size of 120 degrees. The grid search was set up to explore the central region of the crystal surface slab, providing a set of 50 initial configurations for both the solvent and the host molecules. Although aspirin is known to form dimers in its crystal structure, a monomer crystallizing host unit was assumed to interact with the crystal face. This assumption is reasonable linking with the polar solvent environment where dimer pairs of the solute might be easily disrupted. However, there is no experimental evidence, to the best of our knowledge, as to what the detailed solution chemistry of aspirin might be in this solvent, but it has been shown previously for other organic systems that solute dimers are not very stable in polar solvents.43 The configurations built in this way were further optimized using molecular mechanics minimization techniques using the Cerius2 package. The most stable energy after the minimization was taken as the interaction energy of the solvent/ host with the surface. As the energy of interaction of the solvent or the host molecules per unit area (Usolvent/host) is needed to calculate the specific solution surface energy the energies obtained (U′solvent/host) were converted into kcal/m2. To determine an appropriate surface area the reticular areas (SR) of all crystal faces were calculated, where the reticular area is the smallest mesh in the net points (see, e.g., 25), and the minimized energies were taken per unit reticular area, thus:

Usolvent/host )

U′solvent/hostnmolecules SRNA

(14)

where nmolecules is the number of molecules assumed to be binding per reticular area. This model, however, assumes one molecule of each solvent plus one of the host molecular binding per reticular area. Despite this, overall the simulation approach adopted is flexible allowing for an increase of this parameter for larger binding areas and/or smaller molecules. Finally, the interaction energies obtained were scaled with respect to the mole fraction of each component in the system. The final solvent binding energy for each face was then calculated as a sum of the solvent and the host interaction energies and was used to reduce the vacuum specific surface energies derived from the attachment energy calculations (eq 9) using the equation:

- Usolution γhkl ) γvacuum hkl

(15)

The total surface energy of a polyhedral particle was calculated as a function of the polyhedral size expressed with respect to the volume equivalent diameter using the calculated solution specific surface energies for the individual (hkl) faces. From this, the polyhedron having the same surface energy as that calculated for the spherical cluster model based on experimental nucleation studies was found.



Figure 4. Flowchart highlighting route-map and methodology associated with the calculation of the anisotropy of the specific surface energy during nucleation.

Figure 5. Molecular surface packing diagram illustrating the basic approach adopted for surface docking studies as carried out using Cerius2 molecular modeling software and used for surface relaxation and to calculate the interaction energies between the solvents and the host with the host surface crystal growth planes.

4. Results and Discussion 4.1. Determination of Critical Nucleus Size and Specific Surface Energy from Experimental Induction Time Measurements. Measurements of the nucleation kinetics revealed that the specific surface energy for aspirin-38% ethanol solution is 3.99 × 10-3 J/m2 for 15.5 wt% of the solid (see Figure 6). The diameter of the critical nucleus and the number of the molecules expected within it based on the approximation of a spherical nucleus were estimated (eq 4) and given in Table 2. The estimated value for the specific surface energy is in good agreement with other reported specific surface energies for organic compounds. According to Granberg and coworkers,4 organic molecules tend, in general, to show low values of solution specific surface energy. In their work, they calculated the specific surface energy for paracetamol in 20/80% acetone-water mixture as being 2.8 × 10-3 J/m2. Larger values of 4.2 × 10-3 to 8.9 × 10-3 J/m2 from investigation of the nucleation of the polar compound urea in alcohol-water systems have been reported elsewhere,52 con-

Figure 6. Experimental induction time data (3 repeats) associated with batch crystallization of aspirin in mixed ethanol/water solution. Specific surface energy is determined from the slope of the straight line. Table 1. Face Multiplicity and d spacing from BFDH Analysis Together with Attachment and Slice Energies as Calculated Using HABIT Program through Summation of the Dominant Intermolecular Interactions face

multiplicity of the face

d-spacing/ Å

attachment energy/ kcal mol-1

slice energy/ kcal mol-1

{100} {002} {011} {110} {1h11}

2 2 4 4 4

11.37 5.66 5.69 5.70 5.5

-4.49 -9.86 -16.82 -15.46 -16.59

-22.90 -17.53 -10.58 -11.93 -10.81

trasting with much lower values of 0.197 × 10-3 to 0.564 × 10-3 J/m2 for the completely nonpolar alkanes, e.g., as crystallized from diesel fuels. Hence, it can be expected that the solution specific surface energy will be defined with respect to the chemical properties


Hammond et al.

Table 2. Results of the Experimental Induction Time Measurement Showing the Resultant Calculated Critical Nucleation Cluster Size (D) as a Function of Solution Supersaturation (S) Together with the Number of Molecules (N) Per Cluster Calculated Using the Classical Nucleation Theory (eq 4) Based on the Determined Value of the Specific Surface Energy T/°C

σ

D/Å

V/Å3

N

40 41 42 43 45

1.83 1.71 1.60 1.50 1.33

13.07 14.67 16.70 19.30 27.26

1169.28 1655.24 2438.42 3762.04 10609.90

5 8 11 18 50

of the crystallizing material and solvent as well as the interaction between them. Kumar et al.53 reported that dipole-dipole interactions also play a major role in modifying the surface free energy of formation of liquid microclusters at the critical cluster formation stage increasing the critical energy barrier to nucleation. Such effects can be expected to result in greater surface free energies for highly polar molecules compared to the values predicted from classical nucleation theory.53,54 Once the stable nucleus is formed, however, the solution specific surface energy will depend on the solvent affinity for solute molecules at the surface, which will possibly lower the value of the solution specific surface energy reflecting the fact that the unsatisfied intermolecular bonds at the surface would now be compensated via bond formation with solvent molecules. On the other hand, the higher the affinity the solvent molecules have for each other, the less affinity they will have for the surface solute molecules. Therefore, the final value of the solution specific surface energy will inevitably reflect a complex balance between the competing intermolecular interactions between the crystallizing material and its surrounding medium. 4.2. Morphological Prediction. The convergence of the lattice energy summation process as a function of intermolecular distance for aspirin, given in Figure 7a, is associated with a calculated lattice energy of -27.42 kcal/mol and the calculated lattice energy of the minimized crystal structure of -28.14 kcal/mol, which revealed reasonable agreement with the experimental lattice energy of -30.03 kcal/mol55 calculated via eq 11.

Table 3. Calculated Surface Area and the Volumes for Polyhedra with Different Sizes Together with the Surface and Volume Components of the Gibbs Free Energy as Estimated Using Classical Nucleation Theory size Dv/Å 18.40 22.02 25.69 44.04 99.01 110.1

Ds/Å 20.15 24.18 28.21 48.37 108.85 120.94

S/Å 1275 1837 2500 7349 37205 45932

V/Å3 3234 5589 8875 44712 509308 698640

Gs/J 10-19

1.23 × 1.77 × 10-19 2.42 × 10-19 7.09 × 10-19 3.59 × 10-18 4.44 × 10-18

Gv/J -3.92 × 10-20 -6.77 × 10-20 -1.07 × 10-19 -5.42 × 10-19 -6.17 × 10-18 -8.46 × 10-18

Calculated surface attachment energies are given Table 1 with predicted morphology, based on these values given in Figure 7b. The crystal morphology obtained from the morphological simulation was found to be in good agreement with optical micrographs taken of the crystals prepared in this study (Figure 7c) and elsewhere (Figure 7d).56,57 4.3. Estimation of Surface and Volume Free Energy for the Polyhedral Nucleus Shape. The surface and volume free energy relationships derived from classical nucleation theory were used to estimate the size of the critical polyhedral nucleus of aspirin crystals. As can be seen from Figure 1, the critical nucleus diameter corresponds to the maximum point of the function G ) f(D) where D is the size of the nucleus for a given temperature. The calculations were carried out for different temperatures, and the critical sizes of the nonspherical nuclei were predicted. As it is not rational to define the size of an irregular shape with only a single parameter, both the equivalent volume and the equivalent surface diameters were used. The surface area and the volume of the polyhedra of different sizes were calculated and, using eq 1, the surface and volume Gibbs free energies for different sized polyhedra are given in Table 3 and Figure 8a. In this, γ was taken as the experimentally determined value from the nucleation experiments as analyzed using a spherical cluster model with the surface area term S being calculated for the polyhedral model based on the validated crystal morphology. Note the general curve form matches that expected (cf. Figure 1).

Figure 7. Morphological simulation of aspirin: (a) Simulation of crystal lattice energy as function of the interaction distance showing the convergence effect at about 30. (b) Predicted crystal shapes using the attachment energy method. (c) Optical micrograph revealing the typical growth morphology obtained in the crystallization experiment reported here. (d) Published data on growth morphology.57



one particular temperature (313 K) is illustrated in Figure 8b. The calculated values of ∆G versus size were fitted with a function of the kind:

G(D) ) P1D2 + P2D3

Figure 8. (a) Surface and volume free energies as a function of cluster size for polyhedral shaped nucleation cluster at 313 K and at a supersaturation ) 1.83. Both volume and surface area equivalent diameter data are provided; (b) total Gibbs free energy as a function of the nucleation cluster size with the critical size expressed as equivalent volume and surface area diameters which correspond to the extreme in the function G vs size; (c) total Gibbs free energy presented as a function of the nucleation cluster size for different temperatures; inset: critical polyhedral nucleus cluster size, corresponding to the maximum in the Gibbs free energy as function of crystallization temperature increasing in size with the temperature in comparison with the smaller critical spherical nucleus size.

The dependency on the Gibbs free energy calculated using both the volume and surface equivalent diameter models for

(16)

with the first derivative of this function giving the critical size, revealing that the diameter of the nucleus at 313 K estimated with respect to the volume and surface equivalent size is 16.54 and 17.56 Å, respectively. This result shows that the size of a polyhedral nucleus will be larger than a hypothetical spherical nucleus calculated from the measured induction data (e.g., D ) 13.07 Å see Table 2) under the same conditions. Mindful of the expected closeness for the calculated values for the volume and surface equivalent diameters, the volume equivalent model was used hereafter. The increase of the critical nuclei size with temperature is shown in Figure 8c in comparison with the related change of the calculated critical nuclei as determined from the experimental data using the spherical model. In this, the larger polyhedral critical nucleus size compared to that calculated on basis of spherical cluster is expected, reflecting the fact that the work required to form this more complex shape from liquid phase is likely to be larger. 4.4. Modeling Anisotropy of the Solution Specific Surface Energy. The dimensionless surface areas (Si/Stotal) of the various crystal surfaces in the crystal habit of aspirin are given in Table 4 column 2. The specific surface energies, calculated from the attachment energies, given in Table 4 column 3, provide an upper bound on the specific surface energy for this material (i.e., under vacuum conditions). Hence, calculations of the specific surface energy using eq 12 reveal a significantly larger value, 93.4 × 10-3 J/m2, compared to that experimentally derived, 3.99 × 10-3 J/m2. The solvent interaction with the nucleus surfaces calculated from the molecular modeling studies are given in Table 4 columns 4-6 for the three molecular species associated with this crystallization system, i.e., ethanol, water, and aspirin. Calculation of the reticular areas for all the dominant faces of aspirin revealed the predominant habit faces, {100} and {002}, to have small reticular areas of 75.10 and 75.33 Å2 compared to less morphologically important habit faces reflecting their dense crystal surface packing. The modeled interaction energies were scaled in accordance with the relative concentration of the solute and solvent molecular species obtained from the solution solubility as used in the experimental determination of the critical cluster size. This provides a measure of the relative probabilities of an encounter between a given crystal surface with, respectively, an aspirin, water, or ethanol molecule. In this work, the mole ratio for aspirin (nh), ethanol (ns1), and water (ns2) was 0.086, 0.652, and 3.4, respectively. From this, the resulting interaction energies (Usolvent1, Usolvent2, Uhost), due to the effect of these species in the solution for each particular face (hkl), were calculated using the equation:

Usolution ) ns1Usolvent1 + ns2Usolvent2 + nhUhost

(17)

The values obtained from eq 17 were used to amend the interfacial energies calculated in a vacuum (eq 9) to evaluate the specific surface energy for all nucleus faces taking into account the effect of the solvent (Table 4 column 7). The results show that the {100} face has the lowest solution specific surface energy compared to all the other crystal faces. This result is quite rational, reflecting the fact that the surface


Hammond et al.

Table 4. Dimensionless Surface Areas of the Nucleation Cluster Faces as Calculated Based upon a Polyhedral Shapea

face

dimensionless surface area, Si/STOTAL

specific surface energy in a vacuum/Jm-2

{100} {002} {011} {110} {1h11}

0.2316 0.0942 0.0547 0.0218 0.0109

-0.0829 -0.0904 -0.1550 -0.1430 -0.1480

interaction energy/kJ mol-1 ethanol water host -31.51 -83.66 -30.09 -37.29 -37.62

-35.24 -25.91 -42.02 -28.17 -29.67

-80.14 -66.83 -61.77 -71.15 -75.04

specific surface energy; result from the model/Jm-2

scaled specific surface energy from the model/Jm-2

0.00365 0.01090 0.10990 0.10800 0.11310

0.00348 0.01040 0.10490 0.10309 0.10795

a The attachment energy calculation has been carried out with HABIT98 program. Interaction energies have been calculated with Cerius2 after minimizing 50 initial configurations solvent/host-surface generated through grid search methods. The face specific surface energies have been scaled with respect to the factor γmodel/γexperiment.

Figure 9. Molecular surface packing diagram showing the chemistry of {100} face (left), which is rich in H and O atoms from the carbonyl and methyl groups which can be expected to influence the interaction with a solvent medium through the formation of hydrogen bonds and hence the surface energy in comparison with the {002) face (right) which is characterized by having aliphatic H atoms from the phenyl rings, resulting in weak interacting with the polar solvent and hence a higher value for the specific surface energy.

Figure 10. The total surface energy of a polyhedral particle presented for different sizes of the polyhedron. A polyhedral nucleus having the same total surface energy for any temperature as a spherical nucleus will have a smaller size compared to the spherical one.

chemistry for this face is characterized by the intermolecular hydrogen bonds aligned perpendicular to surface plane, hence making it strongly hydrophilic and rich in donor-acceptor interactions.58 This is especially true for crystallization from a polar solvent medium where this particular surface chemistry would cause a decrease in the specific surface energy of this surface due to the effects of additional solvent interactions. The other dominant face {002} exposes mainly the H atoms from the phenyl groups arranged at the surface, which do not interact that strongly with a polar solvent, hence making this face far more hydrophobic and therefore resulting in a higher solution specific surface energy calculated (Figure 9). 4.5. Confrontation of the Model with Experimental Data. The calculated values of the solution specific surface energies for the individual habit faces were then substituted in eq 12 to estimate the total polyhedral specific surface energy (γmodel TOTAL), which was calculated to be 41.9 × 10-3 J/m3 somewhat higher

when compared to that directly obtained from the nucleation experiments (3.99 × 10-3 J/m3). Figure 10 shows a plot of the total solution surface energy calculated as a function of volume equivalent diameter for both the polyhedral shape and that determined from experimental data using a spherical nucleus. From this correlation, it can be seen that the size for the polyhedral particle would be ca. 3 times smaller compared to that based on a spherical cluster; this is assuming that the total surface energies for the two cluster models are equal. This discrepancy in these two sizes reflects the fact that molecular modelling calculations over-predict the total specific surface energy in comparison with the value extracted from the experimental data. When the approximations inherent in the current surface modeling methodology adopted here are kept in mind, these results are in rather pleasing agreement and show significant improvement with respect to other modeling studies. This is not unexpected, reflecting the molecular and crystallographic (system-specific) basis for the modeling approach developed and used here. Hence to obtain representative values of the individual solution specific surface energies the values obtained from the model were scaled via:

γ experiment

γscaled ) γimodel/ i

model γTOTAL

(18)

with results given in Table 4 column 7. However, it is envisaged that further improvements to this new methodology used for the calculation of surface energies might significantly reduce any remaining discrepancy. For example, the specific surface energy obtained from the polyhedral model would be further reduced if the binding of more than one molecule per reticular area (as currently assumed by existing model) was assumed. This would be expected especially to affect the smaller and less dominant habit faces, which have the larger reticular areas. However, simulation of the crystal/solution interface in great detail would require application of more rigorous molecular dynamics


modeling methods,59 contrasting with the methodology developed and described here which is much simpler and more flexible.

Crystal Growth & Design, Vol. 6, No. 6, 2006 1333 Supporting Information Available: Calculation of Cartesian coordinates of the vertices of an irregular polyhedron; calculation of the volume and surface area of an irregular polyhedron; cross product of two vectors lying on a plane (figure of hkl face). This material is available free of charge via the Internet at http://acs.pubs.org.

5. Conclusions In this work, a molecular perspective to modeling solutionmediated nucleation processes has been developed and tested through the evolution of a detailed modeling framework for calculating many of the parameters that would be expected to influence nucleation including solution specific surface energies, surface and volume free energy, critical nuclei cluster size, and solvent binding effects. The model is based on presenting the nucleus as an irregular polyhedron reflecting the final crystal growth morphology, which avoids the unrealistic assumption of a spherical nucleus as applied in classical homogeneous nucleation theory. This model for the crystal nucleus can be expected to be applicable for compounds that do not undergo any polymorphic transformation following nucleation and growth of the crystal, e.g., for a system such as aspirin as studied here. Interestingly, recent studies on L-glutamic acid60 have revealed small clusters post-nucleation to have different polymorphic stabilities. It is worth noting that the methodology developed is independent of the method used to predict the nucleus shape. The latter could be, of course, improved with using a more sophisticated approach for calculating cluster shape, for example, as a function of the solution supersaturation. However, a clear assumption in the approach adopted so far is that there is no significant change in particle shape when going from initially formed crystal nucleus through to the final crystal size. The model has been applied to a test case example based on crystallization of aspirin from mixed aqueous/ethanol solution. The crystal morphology of aspirin has been simulated using attachment energy modeling from which the Cartesian coordinates of the vertices of a polyhedral crystal were derived. The results reveal a prismatic morphology for aspirin crystals consistent with previous data and experimental validation of the shape presented here. The surface areas and the volume of this shape, as estimated based on the polyhedral form, reveal that the anisotropic polyhedral nucleus has larger critical cluster size for nucleation compared to the isotropic spherical cluster model for the same process conditions of temperature and supersaturation. Calculations of surface energies based on attachment energies have been adjusted to allow for the effect of the solvent binding to the habit surface as modeled using the surface docking approach.47,48 The values obtained are consistent with the surface chemistry of the crystal faces but slightly overestimate the total crystal/solvent interfacial energy compared to that for a spherical nucleus as calculated from solution phase experimental nucleation data. In addition, although the model does not take into account the effect of sizedependent particle shape characterization, this overall approach nevertheless provides a useful and potentially important model platform for subsequent future development. Acknowledgment. This work, which forms part of the Ph.D. studies of one of us (K.P.), has been carried out as part of an EPSRC collaborative programme with Dr. Peter Halfpenny at University of Strathclyde in Glasgow, UK, which has been funded by EPSRC Grants GR/R/14491 and GR/R/19328. We also gratefully acknowledge Dr. Liviu Marin, School of Mathematics, University of Leeds, for helpful discussions concerning the surface area and volume calculations. One of us (K.P.) gratefully acknowledges the ORS and Tetley and Lupton scholarship schemes for funding support.

Nomenclature A0 B0 Dcritical Dv Ds dhkl Eatt Elatt Es Esl GS Gv Hsub ∆G ∆Gv ∆Gs J KB M Mi NA ns(h) nmolecules P1,2 p R Rgrowth S Si Ss SR Stotal T V Vcell Usolvent/host U′solvent/host Usolution Z γ γi γimodel γTOTALmodel γexperiment ν σ τ SPGR

preexponential factor preexponential factor, reciprocal of A0 diameter of the critical nucleus, m equivalent volume diameter, m equivalent surface diameter, m d spacing, Å attachment energy, kcal/mol lattice energy, kcal/mol surface energy, J slice energy, kcal/mol surface free energy, J volume free energy, J enthalpy of sublimation, kcal/mol Gibbs free energy change, J volume Gibbs free energy change, J surface Gibbs free energy change, J nucleation rate, number/s m3 Botzman constant 1.3805 × 1023, J/K transformation matrix multiplicity of crystal face Avogadro number 6.023 × 1023 mol-1 mole fraction of solvent (host) number of molecules per reticular area fitting coefficient half of perimeter of a triangle, m gas constant 8.314, J/mol K relative growth rate surface area, m2 surface area of particular crystal face, m2 surface entropy, J/K reticular area, m2 total surface area of particle, m2 temperature, K volume, m3 crystal cell volume, A˙ 3 interaction energy of solvent/host molecule with crystal face, kcal/m2 interaction energy of solvent/host molecule with crystal face, kcal/mol interaction energy of all solution species with crystal face, kcal/m2 number of molecules in unit cell specific surface energy J/m2 specific surface energy of particular crystal face J/m2 specific surface energies for individual faces as predicted from the model, J/m2 total specific surface energy as predicted from the model, J/m2 total specific surface energy as calculated from experimental data, J/m2 molecular volume, m3 supersaturation induction time, min space group

References (1) Volmer, M. J. Phys. Chem. 1926, 119, 277. (2) Kaldis, E.; Scheel, H. Crystal Growth and Materials; North-Holland Publishing Company: New York, 1977; Vol 2. (3) Larson, M.; Garsidel, J. J. Cryst. Growth 1986, 76, 88. (4) Granbergl, R.; Ducreux, C.; Gracin, S.; Rasmuson, A. Chem. Eng. Sci. 2002, 56, 2305. (5) Rajendran, K.; Rajasekeran, R.; Jayaraman, D.; Jayavel, R.; Ramasamy, P. Mater. Chem. Phys. 2003, 81, 50. (6) Roberts, K.; Sherwood, J.; Stewart, A. J. Cryst. Growth 1990, 102, 419.

1334 Crystal Growth & Design, Vol. 6, No. 6, 2006 (7) Hudson, A. AdV. Colloid Interface Sci. 1984, 21, 303. (8) Kashchiev, D. Nucleation - Basic Theory with Applications; Butterworth Heinemann: Oxford, 2000. (9) Oxtoby, D. Acc. Chem. Res. 1998, 31, 91. (10) Horst, J. H.; Jensen, P. J. Surf. Sci. 2005, 574, 77. (11) Horst, J. H.; Kramer, J. M.; Jensen, P. J. Cryst. Growth Des. 2002, 2, 351. (12) Anwar, J.; Boateng, P. J. Am. Chem. Soc. 1998, 120, 9600. (13) Shore, J.; Perchak; D.; Shnidman, Y. J. Chem. Phys. 2000, 113, 6276. (14) Nyvlt, J. J. Cryst. Growth 1968, 3, 377. (15) Mullin, J. Crystallization, 4th ed.; Butterworth Heinemann: London, 2001. (16) Myerson, A. Handbook of Industrial Crystallization; ButterworthHeinemann: New York, 1993. (17) Seville, J. Processing of Particulate Solids; Blackie Academic and Professional: London, 1997. (18) Chenthamarai, S. J. Mater. Sci. Lett. 1996, 15, 1587. (19) Van Gelder, R., Ph.D. Thesis, University of Strathclyde, Glasgow, UK, 1998. (20) Roberts, K.; Walker, E.; Glydesdale, G. The Crystal Habit of Molecular Materials: A Structural Perspective, In Theoretical Aspects and Computer Modeling of the Molecular Solid State; Gavezzotti, A. Ed.; Wiley: New York, 1997. (21) Hartman, P.; Bennema, P. J. Cryst. Growth 1980, 49, 145. (22) Hartman, P. Acta Crystallogr. 1955, 8, 49. (23) Docherty, R.; Glydesdale, G.; Roberts, K.; Bennema, P. J. Phys. D Appl. Phys. 1991, 24, 89. (24) Docherty, R.; Roberts, K. J. Cryst. Growth 1988, 88, 159. (25) Phillips, F. Introduction to Crystallography; Longmans: Bristol, 1955. (26) Roberts, K. Final report to EPSRC, Grant GR/H/40891; University of Strathclyde, Glasgow, 1995. (27) Clydesdale, G.; Docherty, R.; Roberts, K. In Controlled Particle, Droplet and Bubble Formation; Wedlock, D., Ed.; Butterworth Heinemann: Oxford, 1994. (28) Clydesdale, G.; Tefler, G.; Roberts, K.; Grant, D. J. Pharm. Sci. 1997, 86, 135. (29) Kuznezov, V. Surface Energy of Solids; Majesty Stationary Office: London, 1957. (30) Van Der Voort, E.; Hartmann, P. J. Cryst. Growth 1988, 89, 603. (31) Villalba, G.; Landrove, R.; Guaregua, J. Langmuir 1997, 13, 1644. (32) Walton, A.; Whitman, D. J. Chem. Phys. 1963, 40, 2722. (33) Walton, A. J. Chem. Phys. 1963, 39, 3162.

Hammond et al. (34) Zaitseva, N.; Rashkovich, L.; Bogotyreva, S. J. Cryst. Growth 1995, 148, 279. (35) Longuemard, P.; Jbilou, M.; Hermann, A.; Guyot, J. Int. J. Pharm. 1998, 170, 51. (36) Clydesdale, G.; Docherty, R.; Roberts, K. Comput. Phys. Commun. 1991, 64, 311. (37) Momany, F.; Garruthers, L.; McCuire, R. J. Phys. Chem. 1974, 78, 1595. (38) Walker, E.; Roberts, K.; Maginn, S. Langmuir 1998, 14, 5620. (39) George, A.; Harris, K.; Rohl, A.; Gay, D. J. Mater. Chem. 1995, 5, 133. (40) Gay, D.; Rohl, A. J. Chem. Soc.: Faraday Trans. 1995, 91, 925. (41) Gracin, S.; Rasmuson, A. Cryst. Growth Des. 2004, 4. (42) Geertman, R.; Heijden, J. J. Cryst. Growth 1992, 125, 362. (43) Parveen, S.; Davey, R.; Dent, G.; Pritchard, R. ChemComm 2005, 2, 2. (44) Bond, W. Am. Mineral. 1946, 31, 31. (45) Dowty, E. Am. Mineral. 1980, 65, 465. (46) Cerius2: San Diego, CA, USA, and Cambridge, UK, 1995. (47) Lu, J.; Urlich, J. Cryst. Res. Technol. 2003, 38, 63. (48) Roberts, K.; Walker, E. Mol. Cryst. Liq. Cryst. 1996, 279, 233. (49) Engkvist, O.; Price, S.; Stone, A. Phys. Chem. Chem. Phys. 2000, 2, 3017. (50) Hammond, R.; Ma, C.; Roberts, K.; Ghi, P.; Harris, R. J. Phys. Chem. B 2003, 107, 11820. (51) Goodman, J. Chemical Application of Molecular Modelling; The Royal Society of Chemistry: Trowbridge, Wiltshire, 1998. (52) Hendriksen, B.; Grant, D. J. Cryst. Growth 1995, 156, 252. (53) Kumar, F.; Jayaraman, D.; Subramaniam, C.; Ramasamy, P. Chem. Phys. Lett. 1993, 206, 415. (54) Wright, D.; Caldwell; R.; El-Shall, M. Chem. Phys. Lett. 1990, 176, 46. (55) Kirklin, D. J. Chem. Thermodyn. 2000, 32, 701. (56) Meenan, P. ACS Symp. Ser. 1997, 667, 2. (57) Watanabe, A.; Yamaoka, Y.; Takada, K. Chem. Pharm. Bull. 1882, 30, 2958. (58) Denesh, A.; Davies, M.; Hinder, S.; Roberts, C. Anal. Chem. 2000, 72, 3419. (59) Boek, E.; Briels, W.; van Erden, J.; Feil, D. J. Chem. Phys. 1992, 96, 7010. (60) Hammond, R.; Pencheva, K.; Roberts, K. J. Phys. Chem. B 2005, 109, 19550.

CG0505618

Crystal

Recommend Documents