Modeling Crystal Shape of Polar Organic Materials - American

Modeling Crystal Shape of Polar Organic Materials: Applications to Amino Acids Vered

Bisker-Leib†

and Michael F.

Doherty*,‡

Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003, and Department of Chemical Engineering, University of California, Santa Barbara, California 93106 Received June 7, 2002;

CRYSTAL GROWTH & DESIGN 2003 VOL. 3, NO. 2 221-237

Revised Manuscript Received November 1, 2002

ABSTRACT: Models for predicting the shape of organic crystals use geometrical rules and intermolecular interaction energies between the building blocks of the crystal but often ignore the environment in which the crystal grows. The resulting shapes predicted by these models may differ from the actual grown shapes since shape can be influenced by the solvent, impurities, and additional external factors. We present a model to predict the shape of polar organic materials crystallized from solution. The model is based on a Burton, Cabrera, and Frank growth mechanism and accounts for the solute-solvent interactions at the interface. Using this model, we have successfully predicted the shape of amino acids crystallized from solution. Amino acids were chosen as model compounds not only because they are biologically important but also because they are small organic molecules that exhibit various functional groups and have strong electrostatic interactions. Their modeling raised fundamental issues discussed in this article, including the description of charge density, determination of growth unit, selection of a force field, and estimation of crystal-solution interfacial free energy. The shape of R-glycine, when grown from the vapor and from aqueous solution, was predicted and compared with actual growth shapes and with predictions by others. The shape of L-alanine, when grown from aqueous solution, was predicted and compared with the experimental growth shape. Introduction Solution Crystallization in the Pharmaceutical Industry. Over 90% of all pharmaceutical products are formulated in particulate, generally crystalline form.1 The production of particles often involves crystallization from solution, as a separation or purification method and also as a method to generate particles from solution. During crystallization, various physical-chemical characteristics of the substance are defined, including shape and size, chemical purity and stability, bioavailability, solubility, and dissolution rate. Accordingly, solution crystallization has become an important field of research. However, control over the physical form has remained poor, mainly since the basic processes that govern it, nucleation and growth, are not fully understood or controlled.2,3 Different crystal shapes of the same crystalline system might be a result of differences in the rate of growth of the same set of faces or might be a result of the appearance of different faces. The latter case is equivalent to having different moieties of the molecule being exposed on the surface and interacting with the environment. When a substance is crystallized from its vapor, the environment contains only neighboring solute molecules; therefore, the interactions are solute-solute in nature. When the process takes place in solution, the environment contains not only neighboring solute molecules but also solvent molecules and may include impurities and additional objects (vessel walls or dust particles) that can influence the growth. The interactions, therefore, include solute-solute and solutesolvent. * To whom correspondence should be addressed. Tel: (805)893-5309. Fax: (805)893-4731. E-mail: [email protected]. † University of Massachusetts. ‡ University of California.

Different packing configurations of the same molecules produce different unit cells (polymorphs), which is another possible cause for disparity in properties of crystalline products.4,5 Different packing configurations also result in different moieties exposed at the interface; in addition, the density of the face/interface can vary from one polymorph to another. In fact, polymorphs may have such significant differences in their solubilities and rates of dissolutions, that the FDA has increased its regulation on the solid form. The need for a model to predict shape that accounts for external factors has increased over the years. Crystal Engineering, which can be defined as the ability to design a crystalline product with desired characteristics by the successful use of computational models followed by synthesis in the lab, has become an important area of science and technology. To advance this area, there is a need to develop models that are able to simulate the process and predict its outcomes, based on available physical-chemical properties of the solute and the solvent and on process conditions. Model for Predicting the Shape of Amino Acids Crystallized from Solution. There has been a continuous effort to enhance our understanding of the relationship between internal crystal structure and environmental factors (such as temperature, supersaturation level, type of solvent, and impurities) to its external shape. The pioneering work of the mineralogists Bravais and Friedel,6 extended by Donnay and Harker,7 has yielded the BFDH rule. The rule states that the larger the interplanar spacing of a face, correcting for the extinction conditions of the space group, the more prominent the face. The BFDH rule has established that there is a relationship between internal geometrical structure and external morphology.

10.1021/cg025538q CCC: $25.00 © 2003 American Chemical Society Published on Web 12/31/2002

222

Crystal Growth & Design, Vol. 3, No. 2, 2003

The importance of interactions between building blocks of a crystallizing system was first acknowledged by Wells.8 Wells found that the external shape of naphthalene changed from an elongated rod to a flat plate when it was crystallized from cyclohexane and ethyl alcohol, respectively. After a series of experiments for additional crystalline systems (iodoform, pentaerythritol, and anthranilic acid) were conducted, he concluded that “all crystals from a solution in a given solvent, if crystallized sufficiently slowly, do exhibit a standard shape which is characteristic of the particular solutesolvent combination.” Wells related the shape to the direction(s) of the strongest interatomic or intermolecular bonds. Hartman and Perdok further developed the idea of strong bonds9a,b in their Periodic Bond Chain (PBC) theory. Their hypothesis is that periodic repeated chains of strong nonbonding interactions determine the external shape of a crystal. They classified crystal faces according to the number of PBCs that were found in the plane parallel to a face (as two, one, or none for F (flat), S (step), or K (kink) faces, respectively). F faces are the slow-growing faces, followed by the S faces, while K faces grow rapidly. A series of contributions by Hartman and Bennema9c,d provided the Attachment Energy model, which was the first model based on a qualitative growth mechanism, in the form of layers. The energy required to form a layer in the bulk is defined as the slice energy, and the energy released upon an attachment of the slice to a face of a growing crystal is the attachment energy. The basic premise of the model is that faces that interact strongly (those with large negative attachment energies) will grow fast relative to faces that interact weakly (those with less negative attachment energies). Therefore, the faces that will appear on the final form are the ones with small absolute attachment energies. The development of the attachment energy model was an important milestone in the prediction of crystal morphology. In the absence of a solvent, or in systems where the solvent does not interact strongly with the solute, the model provides good predictions.9c,d,10,11 The use of the attachment energy model to predict crystal shape is most appropriate for crystals that are grown from the vapor. When grown from solution, the resulting shapes may resemble the actual grown shapes, but in many instances, they differ.10-14 Dissimilarities are attributed to the fact that the model does not account for external factors, such as solvent or impurities, and it is not based on a detailed growth mechanism. Liu and Bennema,15a-e whose model is based on a detailed Burton, Cabrera, and Frank (BCF) growth mechanism and accounts for the solvent, were the first to model in a systematic way the shape of crystals grown from solution. They introduced the concept of effective growth units, which are solute molecules adsorbed at the surface that are in dynamic equilibrium with the solid molecules at the kinks. Their model has two key parameters: thkl, the free energy associated with the transition of an adsorbed solute molecule to an effective growth unit, and Chkl, the surface scaling factor, which accounts for the solvent effect on the surface. Both of them are face-dependent and need to be calculated using

Bisker-Leib and Doherty

Self-Consistent Field (SCF) lattice model calculations or molecular dynamic and/or Monte Carlo simulations. Winn and Doherty14 introduced a model to predict growth morphology, which is based on detailed growth mechanisms (two-dimensional nucleation, BCF mechanism, and rough growth mechanism) and classical interface properties of the solute and the solvent. Their model was applied successfully to nonpolar systems. Utilizing several principles of this model, Bisker-Leib and Doherty16 developed a method to predict crystal shape of either polar or nonpolar organic materials. We will present its major principles and provide some guidance about additional considerations required for the application of this method to amino acids. Amino Acids. Amino acids are named after the ammonium ion, which is a vital functional group in almost all neurotransmitters and in many other important biological molecules. Some amino acids, for example, glycine and γ-aminobutyric acid (GABA), act as essential neurotransmitters in the central nervous system.17,18 Hence, a better understanding of the behavior of amino acids in solution and in the crystalline state can provide insight into their interactions with their receptors, to their mechanisms of action, and to potential disruptions of these mechanisms.19 R-Amino acids, with one amino group and one carboxylic acid group, crystallize from water as dipolar species known as zwitterions or dipolar ions.20d Thus, a zwitterion is an organic molecule, frequently an amino acid, which has a positively charged segment and a negatively charged segment. The existence of a zwitterion depends on its environment, and there is evidence to suggest that it does not exist in the gas phase but is limited to the solution and solid states.20c,21 The nonzwitterionic form of glycine in the gas phase has been confirmed experimentally by means of microwave spectroscopy.22,23 Zwitterions create hydrogen bonds, in the form of N-H+- - -O--C, which are very strong bonds that play a major role in determining the structures of amino acid crystals. Typical hydrogen bonds of water-amino acid systems are in the following formats: -COO-- - -H3+N-, -COO-- - -H2O, and -NH3+- - -OH2. Because the hydrogen bonds created between zwitterions are very strong, their possible structures are confined to a limited number of configurations, i.e., to only energetically stable configurations.20d,21,24 Although hydrogen bonds fall between chemical bonds and nonbonding interactions, in terms of both the range of interatomic distance and the range of energy of interaction,24 we treat them as intermolecular interactions in our calculations. Addadi et al.25 and Weissbuch et al.26 developed a method, based on controlled inhibition of growth through adsorption of tailor-made additives, for the direct assignment of the absolute configurations of resolved chiral molecules. Their method can be applied to polar and enantiopolar chiral crystals. In the case of polar chiral crystals, polar additives were used to selectively inhibit certain faces of the crystal. The observed morphological differences between crystals grown in the presence and in the absence of additives were used to determine the direction of the molecules with respect to the polar axis. In the case of centrosymmetric enantiopolar crystals, they found that the orientation

Modeling Crystal Shape of Polar Organic Materials

Crystal Growth & Design, Vol. 3, No. 2, 2003 223

of the chiral molecules with respect to the crystal axes could be determined by X-ray diffractometery, and they used this technique to probe the configuration of chiral additives. Amino acids such as R-glycine or (R,S)-serine, which crystallize to give centrosymmetric enantiopolar crystals, were used as model compounds for this method. The absolute configuration of the resolved additives was determined through observation of the morphological changes that they selectively induced on the enantiopolar faces of R-glycine or (R,S)-serine. The packing configurations of such similar crystals as chiral-resolved R-amino acids and of their racemic counterparts were found to be the reason for the difference between their inductions of ice nucleation. Ice nucleation has been exploited for the induction of rainfall in clouds, and it is related to damage to nonconiferous plants. Gavish et al.27 proposed a model in which crystals that have a polar axis induce a freezing point higher by 4-5 °C than the corresponding crystals that lack a polar axis. This demonstrates that a better understanding of crystal structure and shape can lead to prediction of valuable characteristics of organic systems. Lattice Energy and Proton Transfer Energy. Lattice energy, Elatt, is the energy of formation of a crystal from the isolated (gas phase) molecules. The lattice energy of molecular crystals equals the sublimation enthalpy, ∆Hsub, if it is corrected by a factor of RRT, where R is normally taken as 2:28a,29

well.20c-d,22,33 In particular, amino acids and peptides, which crystallize as zwitterions, change to neutral form in the gas phase upon sublimation through intramolecular proton transfer. During sublimation, a proton is transferred from a (NH3+) or (NH2+) group to the (COO-) group. The energy required for this process is termed proton transfer energy, and it can only be estimated since it is not known from experiment (for amino acids and peptides, which are nonzwitterionic in the crystal, another deformation energy of about 7.5 kJ/ mol must be taken into account, due to small changes in the geometry of the carboxylic end, which occur in the transition from hydrogen-bonded molecules in the crystal to free molecules in the gas phase22). Although proton transfer energy was already investigated in the 1950s,33 it is not clear how it should be calculated and what its respective values are. No et al.20c proposed a Born-Haber cycle34 in which the relationship between the lattice energy, sublimation enthalpy, and proton transfer energy can be described as follows:

Elatt ) -∆Hsub - RRT

s ∫0T [C g,x zwi (T) - C zwi (T)] dT = RRT

(1)

Docherty30

tested the relationship in eq 1 by plotting the calculated lattice energies of 80 molecular materials against the experimental lattice energies (derived from the negative sum of the sublimation enthalpy plus 2RT). The average difference between calculated and experimental energies was found to be less than 6%. Enthalpy of sublimation can be obtained directly by calorimetry or indirectly using the Clausius-Clapeyron relation (see nomenclature section for definitions of the various parameters):

dP ∆HsubP ) dT RT 2

(2)

If the sublimation enthalpy is independent of the temperature, this equation can be integrated to yield

ln

( )

(

)

-∆Hsub 1 P 1 ) Pref R T Tref

(3)

Therefore, by measurement of the vapor pressure of a solid (e.g., using the Knudsen effusion technique31), its sublimation enthalpy can be determined and its lattice energy can be calculated (because vapor pressures of solids are very small, their measurement is intrinsically difficult. Characteristic vapor pressures of amino acids are in the order of 0.2 × 10-4 to 2.5 × 10-4 Torr at 455 K31,50). This relationship between the lattice energy and the enthalpy of sublimation holds as long as the geometries of the molecules in the crystal and in the gas phase are equal.22,32 However, if these geometries are different, molecular deformation energy has to be considered as

T s ∫0T C g,x zwi (T) dT ) ∫0 C zwi (T) dT + T ∆Hsub(T) - ∫0 C gneu (T) dT + ∆E ab err,neu + ∆Ept T g ∆E ab err,zwi + ∫0 C zwi (T) dT + ∆Erelax(T) (4)

Elatt(0) +

This relationship can be simplified by using the approximation:

(5)

Rae and Mason35 determined that R ) 2, and it has become a standard rule.13,28a,e,29,45 Additional terms are sufficiently small so that they can be neglected (see full development in No et al.20c), which results in (we have reversed the signs of the sublimation enthalpy and 2RT from the original development of No et al.20c to be consistent with the definitions above):

Elatt(0) ) -∆Hsub(T) - 2RT + ∆Ept

(6)

The same result can be obtained by using the development found in Voogd et al.:22

∆Hsub ) H g - H cr ) 4RT + (g0 - cr 0) + g cr (U gvibr - U cr vibr ) + (U el - U el ) - Elatt (7)

Using the approximation of Rae and Mason:35

w Elatt ) -∆Hsub - 2RT + (U gel - U cr el )

(8)

∆Ept ) U gel - U cr el

(9)

Proton transfer energy, ∆Ept, according to Voogd, is the difference between the ab initio molecular energies of the gas (neutral form) and the crystal (zwitterionic form). Voogd et al.22 and No et al.20c carried out extensive ab initio calculations of electrostatic energies and proton transfer energies of several amino acids. Voogd concluded that zwitterionic crystals have larger electrostatic energies than nonzwitterionic crystals; however, the difference is roughly canceled by the proton transfer

224



Table 1. Comparison of Calculated Proton Transfer Energy and Lattice Energy of r-Glycine by Using Different Force Fields and Charge Distributions (the Atomic Numbering Is Given in Figure 1) Voogd22

No20b,c

Kwon20d

Boek12

ab initio 6-31G**

ab initio 6-31G**

ab initio 6-31G** calcd net atomic charges

method

ab initio DZP

point charge method

Mulliken population

Mulliken population

PDa

PEOE

C1 C2 (CR) N3 O4 O5 H6 H7 H8 H9 H10 ∆Hsub (kcal/mol) (measured) ∆Ept (kcal/mol) Eelec (kcal/mol) Elatt (kcal/mol)

0.493 -0.204 -0.467 -0.632 -0.572 0.166 0.170 0.324 0.367 0.355 not reported

0.753 -0.187 -0.624 -0.710 -0.710 0.183 0.183 0.375 0.375 0.375 32.650

0.819 0.065 -0.404 -0.763 -0.763 0.051 0.051 0.314 0.314 0.314 32.650

0.869 0.195 -0.531 -0.785 -0.785 0.012 0.012 0.337 0.337 0.337

-37.28b -49.12 not reported

-34.32c -49.54 -68.0d

-34.32 -40.31 -68.0

-43.1 -64.8

0.1845 0.0394 0.0702 -0.4903 -0.4698 0.0560 0.0531 0.1794 0.1893 0.1879 34.7

-35.85e -73.13

a PD ) potential-derived net atomic charges. b Voodg defines ∆E according to eq 9. c No defines ∆E as the difference between ab pt pt initio gaseous states of the zwitterionic and neutral forms. d No calculated the lattice energy according to eq 6 at room temperature. e Boek calculated this value and then added proton transfer energy from Voogd et al.22

energy. According to Voogd’s development, the contribution to the lattice energy is the difference between the electrostatic energy and the proton transfer energy. This difference approximately equals the electrostatic energies of nonzwitterions (Table 4 in Voogd et al.22) and, therefore, leads to comparable enthalpies of sublimation. This may account for the similarity in sublimation enthalpies of zwitterions and nonzwitterions and to the observation that some amino acids and peptides crystallize as zwitterions and others as nonzwitterions. Voogd et al.22 calculated electrostatic lattice energies and proton transfer energies of various amino acids using quantum mechanical methods with four ab initio basis sets in order of increasing quality. Electrostatic interactions were calculated using point charges obtained from gross Mulliken populations. The proton transfer energies calculated were (surprisingly) high, in the order of neutral sublimation enthalpies, and showed high dependency on the basis set selected. For instance, proton transfer energies for R-glycine were estimated at -39.9, -32.7, and -37.3 kcal/mol for the basis sets SV, DZ, and DZP. Similar differences were found for additional amino acids. The calculation of electrostatic lattice energies revealed a similar trend, as electrostatic energies of R-glycine were estimated at -49.3, -52.1, and -49.1 kcal/mol according to the basis set selected. The DZP basis set (last one quoted) is generally regarded as the most accurate set and leads to the most reliable value73 for the calculated electrostatic energy of -49.1 kcal/mol. No et al.20c calculated electrostatic energies and proton transfer energies of neutral and zwitterionic amino acids using an ab initio SCF method with the 6-31G** basis set. His results for glycine are reported in Table 1 and compare well with Voogd’s results using an ab initio method with the DZP basis set. Our model deals with growth from solution at the solvent-solute interface, where the molecules are already in the zwitterionic state. The zwitterionic state is usually energetically favored. For instance, at room temperature, the zwitterionic state of glycine is favored

over the neutral state by a free energy and enthalpy of 7 and 10 kcal/mol, respectively.21,36 Because the molecules in the bulk and at the interface are already in a zwitterionic state, we hypothesize that their interactions can be represented by constructing a crystal model from the zwitterionic crystallographic data, applying an appropriate force field and a charge description. Consequently, there is no need to account for proton transfer, which happens in the gas phase. However, for growth from the vapor, this energy should be taken into account. We are aware of one attempt12 to incorporate proton transfer energy into models for determination of shape of zwitterions. Description of the Intermolecular Interactions. A reliable description of the intermolecular interactions between the constituent units of a molecular crystal is a prerequisite for lattice energy calculation and a basic requirement for most models that can be used to predict crystal shape, packing configurations, protein conformations, etc. Over the last four decades, there has been a considerable effort to describe in a precise manner the interactions between molecules in organic crystals.20,22,28,32,37-42 The dominant intermolecular forces can be described as Van der Waals (VdW), electrostatic, hydrogen-bonding, and specific directional interactions. VdW interactions have attractive and repulsive parts. The attractive part (London dispersion) can be considered to result from the small transient dipoles that molecules can induce in each other. Attractive forces show a common dependence of the distance between the atoms to the sixth power. Repulsive forces occur as most stable molecules have closed electron shells and cannot accept other electrons without violating the Pauli exclusion principle. Repulsive forces decrease with the ninth to the twelve power of the interatomic distance (sometimes described by an exponential dependency), and they dominate at very short distances. The sum of the attraction and repulsion energies (VdW interactions) is generally regarded as nonbonded energy.



Electrostatic interactions can be described in the form of a multipolar expansion. The decrease of the interaction energy with the order of the multipole allows neglecting terms beyond quadrupoles. The calculation of multipoles is not straightforward, even for small molecules; therefore, simplified models have gained popularity. In the monopolar model, the total energy is calculated as the sum of all pairwise products of the atomic partial charges divided by the distance between the centers of the atoms. Hydrogen bonds can be qualitatively defined as electrostatic interactions. Practically, they can be calculated by dividing the hydrogen bond interaction into two terms: the “usual” nonbonded energy and an electrostatic term that takes into account the partial charges on the donor and acceptor groups.11,29,43 The orientation of nonpolar organic molecules in a crystal results from a balance between VdW attraction and repulsion forces. VdW interactions also exist between polar organic molecules, but electrostatic interactions and hydrogen bonds play key roles. Specifically, hydrogen bonds greatly impact the packing configuration of molecules, as observed packing is almost inevitably that allowing the maximum number of such bonds to be made.29 A comprehensive review on the role of hydrogen bonding in molecular assemblies is provided by Bernstein, Etter, and Leiserowitz.24 Combining the attractive and repulsive forces with a simple Coulombic potential to account for the electrostatic interactions, one can obtain a description of the intermolecular interactions. The two expressions most commonly used to represent combined VdW and Coulombic interactions are the Lennard-Jones 6-12:

V LJ )

qiqj -A B + 12 + 6 r r r

(10)

and Buckingham 6-exp functions:

V X6 )

qiqj -A + B exp(-Cr) + 6 r r

(11)

where A, B, and C are empirical atom-atom parameters, qi and qj are the partial charges placed on atoms i and j, and r is the interatomic distance. A, B, and C can be estimated from adjusting the chosen potential to observed properties, including crystal structures, heats of sublimation, and hardness measurements (ref 28e,f). An additional term may be added to the above descriptions to account for hydrogen bonds explicitly. Kitaigorodskii37a-e and Williams38a-e were the first to derive parameters for describing atom-atom interaction and to use them for crystal structure packing analysis. Their approach, which is called the atom-atom approximation, assumes that the interaction between two molecules consists of the sum of all of the nonbonded interactions between the constituent atoms. If there are n atoms in the central molecule and n′ atoms in each of the N surrounding molecules, then one gets N

Elatt ) 1/2

n

n′

∑ ∑ ∑ vkij k)1 i)1 j)1

(12)

where vkij is the interaction energy between atom i in

the central molecule and atom j in the kth surrounding molecule. Because the interactions are between two bodies, the factor of 1/2 arises to avoid double counting. William’s38d,e set of VdW parameters, which is based on neutron diffraction data and available sublimation enthalpies, is considered to be highly accurate, and it is still used with minor changes in recently developed force fields, such as Dreiding.41 Calculation of Net Atomic Charges. In molecular crystals that have polar groups, electrostatic interactions predominate at intermediate and larger distances. Moreover, electrostatic interaction is sensitive to the orientation of the molecules and, therefore, plays a major role in determining the packing structure and the shape of crystals.28a,e,44 Because of the importance of electrostatic interactions in determining the shape, a precise description of molecular electron density is a prerequisite for a successful model of the system. Compactness is another desired feature in the description of electron density, to allow a simple simulation using computational methods.13,20a However, “although chemists have an intuitive feeling for the qualitative nature of charge distributions the assignment of quantitative values is met with great difficulties.”47 Molecular electron density is a well-defined quantum mechanical property, and the electrostatic potential can be calculated from it.35,48b Although generally agreed that the electron density can be approximated by a superposition of atom-centered point charges, it is not clear how such net atomic charges should be defined. Numerous quantum mechanical methods have been devised as solutions to this problem. They can be classified, in order of decreasing accuracy, as methods involving distributed multipolar expansions, multipoint charge distribution, and point net atomic charge distributions. The best approximation is achieved by representing the molecular electron density by a set of atomic multipoles (charges, dipoles, quadrupoles, etc.), since it accounts in a very precise way for the anisotropic nature of molecular charge distribution.20a In the multipoint distribution method, monopoles are positioned within the molecular scaffold at places other than atomic centers. However, the simplest (thus common) method for describing the distribution of charge in a molecule is to use atom-centered point charges. Improved quantum mechanical methods include calculating ab initio molecular electrostatic potential and fitting a set of atomic monopoles that best reproduce such a potential. The potential-derived (PD) method, developed by Cox and Williams, is an example of such a method. Charges are determined by least-squares minimization of the differences between the quantum mechanical electrostatic potential, calculated ab initio, and the potential generated from the atomic charges, on some points on a grid of regularly spaced points surrounding a molecule. Charges derived by this method usually reproduce the permanent electric moments of a given molecule calculated from the wave function better than a population scheme, but the method is applicable only to small molecules, for which ab initio calculations are feasible.20a-c,42 Alternatively, the electron density can be determined by experimental methods, including X-ray diffraction

226


experiments at very low temperature. These experiments provide density maps, from which atomic net charges, dipole moments, and quadrupole moments can be derived and used for the parameters in the electrostatic term of the potential energy.46 Gasteiger and Marsili47 developed a different method to calculate partial charges, Partial Equalization of Orbital Electronegativity (PEOE), based on the concept of orbital electronegativity introduced by Hinze and Jaffe.49 According to Gasteiger and Marsili,47 orbital electronegativity can be described in terms of ionization potentials and electron affinities. Their iterative procedure transfers the charge along the covalent framework of the molecule until calculated point charges reproduce electrical multipole moments and ab initio electrostatic potential. Upon transfer, an electrostatic field is generated, which makes further charge transfer more difficult; therefore, charge transfer is blocked before total equalization of electronegativity is reached. This method was modified by No et al.20a,b In the modified method, electronegativity is calculated using parameters based on gas phase experimental dipole and quadrupole moments. Additionally, in the original formalism, a single damping factor, f, was introduced, to simulate the effect of the electrostatic field on the transfer of charge. In the modified method, so-called MPEOE, this damping factor is replaced with a set of five damping factors, which were tailored to fit the different orbitals. An attractive feature of both methods (PEOE and MPEOE) is that they involve direct calculation of the charge distribution for a given molecule using available atomic and orbital characteristics. Unfortunately, the point charges calculated by the above methods are not transferable from one molecule to another, even for atoms located in similar chemical environments; therefore, a separate calculation needs to be performed for every molecule. Because the description of charge distribution was found to have a great influence on the resulting calculated solid state intermolecular interactions and hence the shapes, we used partial charges calculated by the PD method. We report a comparison between calculations of lattice energy, proton transfer, and electrostatic energy of R-glycine (Table 1) and L-alanine (Table 2) based on different force fields and charge distributions. Determination of the Growth Unit. In modeling of crystalline systems, a specification needs to be made about the nature of the crystal precursor, the basic growth unit.48a-c,51 A growth unit can be described as the smallest cluster (a molecule, dimer, trimer, etc.) that is created in solution during nucleation and docks onto crystal faces during crystallization. For modeling purposes, it is considered a rigid body, interacting with the surrounding clusters by intermolecular interactions, which ultimately influence the resulting shape. Thus, a determination of a growth unit has further implications; among them, it includes the classification of certain interactions as nonbonding or bonding interactions. The concept of a growth unit should not be confused with a different concept, an asymmetric unit,48a which is the smallest structural unit that upon translation in space according to the symmetry operations will reproduce the periodic structure of the crystal. Even though the structure of the growth unit is

Bisker-Leib and Doherty Table 2. Comparison of Calculated Proton Transfer Energy and Lattice Energy of L-Alanine by Using Different Force Fields and Charge Distributions (the Atomic Numbering Is Given in Figure 2) Voogd22

No20b,c

No20b

method

minimal basis

ab initio 6-31G**

ab initio 6-31G**

point charge method

Mulliken population

Mulliken population

PDa

modified PEOE

0.783 -0.101 -0.351 -0.610 -0.721 -0.721 0.368 0.368 0.368 0.179 0.146 0.146 0.146 33.050

0.883 0.245 -0.359 -1.199 -0.772 -0.772 0.516 0.516 0.516 0.071 0.118 0.118 0.118 33.050

0.877 0.149 -0.014 -0.412 -0.783 -0.783 0.303 0.303 0.303 0.020 0.012 0.012 0.012 33.050

C1 0.486 C2 (CR) 0.015 β C3 (C ) -0.372 -0.805 N4 O5 -0.654 O6 -0.620 H7 0.423 H8 0.455 H9 0.445 H10R 0.176 β H11 0.181 H12β 0.145 H13β 0.125 ∆Hsub (kcal/mol) not reported (measured) ∆Ept (kcal/mol) -37.76b Eelec (kcal/mol) -56.17 Elatt (kcal/mol) not reported

-31.34c -47.12 -65.5d

-31.34 -42.70 -65.5

-31.34 -43.79 -65.5

a PD ) potential-derived net atomic charges. b Voogd defines ∆Ept according to eq 9. c No defines ∆Ept as the difference between ab initio gaseous states of the zwitterionic and neutral forms. d No calculated the lattice energy according to eq 6 at room temperature.

important for modeling purposes, some studies tend to overlook the importance of its determination and others specify a certain growth unit based on intuition. There are systems, e.g., benzamide10a,b and propionic acid,28e where the growth unit can be identified as a dimer, based on the orientation of the molecules in the solid state. In other systems, e.g., R-glycine10a,b and benzoic acid,52 it is not clear whether the growth unit is a monomer or a dimer. Berkovitch-Yellin,10a who crystallized R-glycine from aqueous solution and from the gas, observed the absence of the large {010} face from the solution-grown crystals that appeared on the vaporgrown crystals. She hypothesized the existence of a precursor in the form of a cyclic hydrogen-bonded dimer. Support of her hypothesis was found in diffusion coefficient measurements of R-glycine by Myerson et al.,53,54 which yielded an average cluster size of ∼1.8 molecules in supersaturated solution. More recently, grazing incidence X-ray diffraction (GID) experiments by Gidalevitz et al.55 determined that for R-glycine grown out of water, the essential building block is a centrosymmetric dimer. The physical orientation of the molecules in the crystal can be used as an indication of the nature of the growth unit. For example, in the case of propionic acid, the molecules are packed in dimer rings motif and in the case of benzamide, hydrogen-bonded cyclic dimers are interlinked along the b-axis to form hydrogenbonded ribbons, which are stacked along the a-axis primarily due to VdW interactions.10b However, the existence of hydrogen-bonded pairs in the solid state does not rule out the possibility of a monomer growth unit, as investigated by Grimbergen and Bennema52 for the case of benzoic acid. Benzoic acid forms hydrogenbonded pairs in the solid state. Grimbergen and Bennema predicted the shape of benzoic acid using two


growth units, a dimer and a monomer. Although the two morphologically important faces, {200} and {001}, were predicted by both types of growth units, neither of the predictions yielded the complete experimental shapes. We have found that for R-glycine, different growth units result in the prediction of different faces, as was observed for benzoic acid.52 Therefore, it is essential to establish an appropriate growth unit. One way of predicting the growth unit is to study solute precursors in solution, as proposed by Gavezzotti et al.48c They used molecular dynamics to investigate the possible configurations of two molecules of tetrolic acid in a solvent box containing 226 carbon tetrachloride solvent molecules. They showed that the cyclic hydrogenbonded dimer was the most persistent configuration. Recent studies of solute precursors in solution51,56 have made significant progress toward the prediction of crystal structures. However, it is possible that small clusters assume different structures when passing from solution to the faces of growing crystals.48c We suggest selection of the monomer as the initial trial growth unit, since it is a basic structural unit that is likely to be the growth unit of nonpolar solutes in solution. In the next step, other growth units might be evaluated based on the solid state structure and the bonding nature of the crystal. Hydrogen-bonded systems are more likely to have growth units in the form of clusters due to the strength of the bonds. Determination of an Appropriate Force Field. Modeling of a crystal involves the selection of an appropriate force field and a description of charge distribution in the molecule. Several force fields (refs 28a-g) include a set of point charges that can be directly assigned to a molecule. These charges are calculated by fitting observed properties to those calculated from the intermolecular forces, based on the assumption of additivity and transferability of atom-atom interactions. Other force fields (Momany et al.,39a-e Scheraga et al.,40 and Mayo et al.41) do not provide point charges, and a separate calculation of charges is required. The selection of an appropriate force field is challenging for additional reasons. Most force fields are applicable only for a limited set of compounds for which they were created and may not provide reliable results for compounds outside the initial set or outside their defined group (e.g., Hagler Lifson Dauber28e for carboxylic acids, amides, and hydrogen-bonded systems). Some force fields (e.g., Dreiding41), which were created with broad applicability, tend to be too general, in the sense that their results are of moderate precision for larger groups of compounds. Additionally, small differences in interatomic potential may lead to significant changes in the calculated properties of compounds.32 A force field can be used for a certain system once several requirements are satisfied. The force field should include terms for all of the functional groups and bond types in the system, the interactions should converge within a reasonable limiting radius, and the lattice energy should approximately equal the sublimation enthalpy (possibly corrected for zwitterionic and other effects). The approach taken by Momany et al. was to develop a potential energy function for neutral amino acids first39d,e and to extend its applicability to amino acid


zwitterions.39e In their work, the proton transfer contribution was ignored when deriving the potentials. The major drawback in this approach, as was pointed out by Kwon et al.,20d is that calculated binding energies (that correspond to lattice energies) resulted in values that were too small. This is mainly due to underestimation of the electrostatic interactions between amino acid molecules in the solid state. The first force field used is by Nemethy, Pottle, and Scheraga (1983),40 which was developed for amino acids. It is an improved version of a force field by Momany et al.39d,e and was created to include experimental information that became available and a larger set of prototype compounds. These force fields are contributions by a group from Cornell University, led by Prof. Scheraga, and are applicable for crystals of hydrocarbons, carboxylic acids, amines, amides,39a-c and amino acids.39d,e The intermolecular potential energy is taken to be the sum of pairwise interatomic interactions, using the atom-atom approximation. It is partitioned into four contributions: VdW attraction, repulsive nonbonded contribution, Coulombic (electrostatic) contribution, and hydrogen bond contribution. The VdW attractive term is obtained by the Slater-Kirkwood method.43 The repulsive term is calculated by minimizing the lattice binding energy with respect to lattice constants of several different crystalline systems. For the electrostatic term, partial charges were calculated by the molecular orbital CNDO/2 (ON) method.39a-c The hydrogen bond term has a 10-12 form. The hydrogen bond attractive term is taken from previous CNDO/2 computations, and the repulsive term was obtained in the same way as the repulsive term of the dispersive interaction. The ensuing intermolecular potential is a combination of calculated quantities with experimentally adjusted quantities, and it is self-consistent with respect to all contributions to the total potential energy. The second force field used is Dreiding, developed by Mayo, Olafson, and Goddard.40 Dreiding was developed to enable predictions for structures where there are little or no experimental data and for novel combinations of elements. In this force field, atoms are treated according to their type, i.e., hybridization state, implicit number of hydrogens, formal oxidation state, etc. The nonbonded potential energy is expressed as a superposition of VdW, electrostatic, and hydrogen bond terms. For the VdW term, an exponential-6 form is selected. Charges were calculated using the Gasteiger and Marsili47 method (described previously), or alternatively, ignored. However, the inventors recommend the use of PD charges from ab initio calculations, whenever it is possible (i.e., small molecules). Model to Predict the Shape of Polar Organic Crystals. Several growth mechanisms for crystals appear in the literature, but few of them have gained broad acceptability. Evidence suggests that crystal faces grow by a screw dislocation mechanism, by a twodimensional nucleation mechanism, or by rough growth. It is also known that different faces of a crystal may grow by different mechanisms, according to the solutesolvent interactions at the interface (surface free energy) and the level of supersaturation. At low supersaturation levels, or large surface free energies, the screw dislocation mechanism is normally operative. The original

228



theory, developed by BCF in 1951,57,58 proposed that screw dislocations, which exist on crystal faces at low supersaturation levels, provide an infinite source of steps onto which oncoming particles can be incorporated. According to this theory, growth occurs by flow of steps across the surface, which forms a spiral. Spirals were observed on faces of crystals,58-60 but more important, the BCF model circumvents the activation energy barrier for two-dimensional nucleation and thus explains experimentally observed growth rates. At moderate levels of supersaturation, the twodimensional nucleation mechanism may apply. Above a critical level of supersaturation, the face is roughened and growth proceeds at a high rate. The BCF expression for the rate of growth normal to a surface is

vhkl hhkl Rhkl ) yhkl

(13)

BFDH model, there is an inverse relationship between the interplanar spacing and the growth rate; thus, halving the interplanar spacing will lead to a doubling of the growth rate. When the BCF model is applied (eqs 13 and 16), the rate of growth in a certain direction is halved when the interplanar spacing is halved. When polar systems are involved, we propose a method that accounts for solute and solvent interactions in the following way:16 the interfacial free energy of two immiscible phases that interact mainly by dispersive interactions (i.e., do not create polar interactions or hydrogen bonds) can be approximated by the sum of the cohesion energies of the single phases, γl, γs, minus the adhesion energy, where the adhesion energy is approximated by the geometric mean rule:62

γls ) γl + γs - WA ) γl + γs - 2(γdl γs)0.5

(17)

where vhkl is the lateral step velocity, hhkl is the step height, which can be approximated by dhkl for monolayer height, and yhkl is the distance between steps. Because growth occurs in kinks, the lateral step velocity depends mainly on the density of kink sites. The average distance between the kinks on a step coinciding with a PBC is

Note that in the calculation of the adhesion energy only the dispersion force component of the surface free energy of the liquid, γdl , is included. For polar systems, the work of adhesion resulting from hydrogen bonds can be approximated, according to Fowkes,63 by the hydrogen bond component of the surface free energy of the solvent, γhl . The interfacial surface free energy of a polar-polar (solute-solvent) system is thus approximated by

X 0hkl ) ahkl[1 + 0.5 exp(φkink hkl /RT )]

γls ) γl + γs - WA ≈ γl + γs - Kγhl

(14)

where φkink hkl is the energy required to form a kink per mole and ahkl is the intermolecular distance in the direction of the step. The classical ideas of Kaichev61 and Chernov58 were incorporated in a model proposed by Winn and Doherty14 to predict the shape of organic materials when crystallized from solution. They propose that when the BCF mechanism applies, the lateral velocity of a face can be estimated from the density of kink sites (ahkl / X 0hkl ) multiplied by a phkl , the distance the step is propagated by adding a monolayer to it: -1 vhkl ∝ aphkl[1 + 0.5 exp(φkink hkl /RT)]

(15)

The rate of growth of a face is therefore

Rhkl ∝

dhkl p -1 a [1 + 0.5 exp(φkink hkl /RT )] yhkl hkl

(16)

The calculation of yhkl, the distance between steps, and of φkink hkl , kink free energy, is provided in Winn and Doherty.14 The interplanar spacing, dhkl, is calculated by taking into account the type of Bravais lattice system (cubic, hexagonal, rhombohedral, etc.) and the unit cell dimensions. Additionally, it is important to include the extinction conditions of the space group as stated in the International Tables of Crystallography,65 since the same reduction in interplanar spacing due to the existence of screw axes and glide planes that causes morphological extinction (as was discovered by Donnay and Harker) also leads to the extinction of X-ray diffraction patterns. When the extinction conditions affect a certain direction, the interplanar spacing in this direction is halved, and the corresponding growth rate is significantly changed. For instance, according to the

(18)

where K is a constant, normally close to unity. In this work, K ) 1. Procedure for the Application of the Model. The following is a step-by-step procedure, which outlines how to create a model of the crystalline solute, calculate the interactions in the solid state and the surface free energy at the interface with the solvent, and determine the relative growth rates of the crystal faces according to eq 16. (1) Create a lattice model of a crystalline solute. (1.1) Use crystallography data, preferably from neutron diffraction studies, to create the crystalline skeleton. (1.2) For nonpolar solute, use the charges provided by the selected force field or ignore charges when unavailable. For polar molecules, calculate net atomic charges. When possible, use a reliable charge density description from quantum mechanical ab initio calculations. If not possible (due to size limitations or others), use readily available molecular orbital methods to calculate partial charges (e.g., Gasteiger and Marsili47 method). (1.3) Select an appropriate force field. (1.4) Perform calculation of intermolecular solid state interactions, using the monomer as the growth unit. The calculation can be carried out using programs such as Habit45 or Cerius2.66 (1.5) Determine which bond chains should be included in the shape calculations. The first nearest neighbors rule usually applies for nonpolar molecules. For polar molecules, select the shorter polar interactions, which tend to be the largest ones on the list. (2) Apply the BFDH rule to the solute system using cell dimensions, space group symmetry, and extinction conditions as dictated by the space group. The calculation provides an extensive list of all of the likely forms (∼15-50), ranked according to their interplanar spacing. This list can be rapidly reduced by removing all of



the high index faces (those that contain indices >4), since high index faces are unlikely to appear in the final shape. The reduced list (∼ 10-20 forms) is used as a starting point for the model selection process of the actual growth faces. (3) Step 1.4 yields a detailed description of the interactions and indicates which bond chains are parallel to each of the likely faces that were determined in step 2. Analyze the interactions using a graphical interface and the calculated intermolecular interactions. When two molecules (or more) interact exclusively in a strong manner (e.g., hydrogen bonds), select them as the growth unit and rebuild the crystal model with the dimer as the growth unit. Repeat steps 1 and 2. (4) Refine the list of likely faces by removing faces that are not parallel to any bond chain (a semi-PBC analysis). (5) Calculate γkink hkl for each of these bonds using eq 17 for nonpolar interactions and eq 18 for polar interactions. (6) When the growth mechanism is known from experiments, use it for the relative growth rate calculations; otherwise, use the criteria provided by Winn and Doherty.14 (7) Calculate relative growth rates according to eq 16 by using the methodology in Winn and Doherty.14 r-Glycine R-Glycine (H2NCH2COOH) belongs to the R-amino acids (H2NCHRCO2H), which mostly crystallize in the zwitterionic form.24 Because of its small size and simple structure, it can be used as a model compound for study of peptides and proteins, and indeed, it has been extensively investigated.10,12,20a,b,21-26,39e,55,64,67,68 However, only a few of these studies deal with the shape of R-glycine.10,12,67,68 R-Glycine zwitterion (H3N+CH2COO-) crystallizes in the monoclinic space group P21/n.64 Unit cell dimensions of the stable R polymorph of glycine are a ) 5.1054 Å, b ) 11.9688 Å, c ) 5.4645 Å, and β ) 111.697°. It contains four molecules centered at M(1) x, y, z; M(2) x + 1/2, -y + 1/2, z + 1/2; M(3) -x, -y, -z; and M(4) -x + 1/2, y + 1/2, -z + 1/2, as can be seen in Figure 1. The extinction conditions for this space group are65 h00, h ) 2n; 0k0, k ) 2n; 00l, l ) 2n; and h0l, h + l ) 2n. The calculated center of mass for molecule M(1) is at xCM ) 0.1259, yCM ) 0.1179, and zCM ) -0.017. Different values for experimental enthalpy of sublimation are reported in the literature, e.g., ∆Hsub ) 34.7 kcal/mol,12 ∆Hsub ) 32.6 kcal/mol,50 and ∆Hsub ) 31.2 kcal/mol.33 Our calculated lattice energies with the Scheraga force field (Table 3) compare well with these values, when the relationship in eq 1 is considered, while the calculations with the Dreiding force field compare only moderately. Lattice energy calculations showed that the energy converged within r ) 30 Å, which was used to determine N in eq 12. Observation of the bonding situation around R-glycine molecules reveals that the molecules are linked by two relatively short, nearly linear, N-H‚‚‚O hydrogen bonds to form layers parallel to the ac plane. The layers are connected, two by two, by weaker N-H‚‚‚O bonds, forming antiparallel double layers, hydrogen bond bilayers, where hydrophilic groups point to the middle of

Figure 1. R-Glycine crystallizes in the monoclinic space group P 21/n. A unit cell with one of the zwitterion molecules highlighted is shown.

the layer and hydrophobic groups are exposed,24,25 as can be seen in Figure 3. There are no hydrogen bonds in the form of N-H‚‚‚O between neighboring bilayers, but the bilayers are connected by C-H‚‚‚O hydrogen bonds across a glide plane. Donnay-Harker Model. Sublimation of R-glycine, in experiments reported by Berkovitch-Yellin,10a yields crystals with hexagonal plate morphology. The {010} faces are reported as the most prominent, followed by {110}, and the top faces reported are {001} and {101 h }. Shape prediction according to the BFDH model provides a bulky shape, which lacks similarity to actual grown shapes (see Figure 4). Attachment Energy Model. Calculated attachment energies are reported in Table 3 and compared with the results of Berkovitch-Yellin10a and Boek et al.12 The energies were calculated using the program Cerius2 66 with the selected force fields of Scheraga40 and Dreiding.41 Partial charges by No et al.20b,c and by Voogd et al.,22 as reported in Table 1, were used to describe the charge density. All shapes presented in Figure 5 were calculated by the attachment energy model; however, there are noticeable differences in the faces that appear and their relative importance. These differences are due to the various force fields and charge distributions used in the calculations, and they demonstrate the importance of the selection of a force field and the description of charge density. Note the similarity between the shape predicted according to attachment energies of Berkovitch-Yellin and the shape predicted by the BFDH model, and how the other shapes increasingly diverge from the BFDH shape. In particular, face {1 h 01}, which is reported to

230



Table 3. Calculated Attachment Energies for Selected Faces of r-Glycine (kcal/mol) Berkovitch-Yellin10a face

dhkl

charge density force field 020 011 110 101 h 111 h 021 h 120 Elatt

5.98 4.67 4.41 4.36 4.10 3.87 3.72

Boek12

this work

model 1

model 2

model 3

multipolar expansion Hagler25a -6.47 -12.98 -9.02 -11.15

Mulliken Voogd22 Buckingham -13.60 -30.32 -27.26 -44.28

Berkovitch-Yellin Hagler25a -6.52 -13.17 -10.10 -14.47

ab initio 6-31G** Hagler25a -7.19 -24.40 -16.11 -29.24

not reported

-30.12 -73.01

-11.97 -26.42

-19.92 -53.97

Figure 2. L-Alanine crystallizes in the orthorhombic space group P 212121. A unit cell with one of the zwitterion molecules highlighted is shown.

be an important face in Berkovitch-Yellin and also appears in the prediction by Boek et al.,12 does not appear in our predictions. Our predictions are of a {020} plate bounded by {110} and {011} faces but a rectangular rather than a hexagonal shape. Because the sublimation energy reported is ∆Hsub ≈ 33 kcal/mol, the lattice energy is expected to be Elatt ≈ -34 kcal/mol. Boek’s models, as well as a model based on the Dreiding force field with charges by Voogd, do not yield lattice energies in the expected range. Therefore, it is believed that these models do not describe R-glycine in a way that is compatible with experimental observations. The models that are based on the Scheraga force field and two different sets of partial charges yielded the expected lattice energies. Also, attachment energies calculated by these models are similar, which may imply that the selection of charge distribution is not signifi-

PD No20b Scheraga40 -2.95 -14.31 -11.63 -29.87 -25.11 -14.60 -12.24 -36.43

Mulliken Voogd22 Scheraga40 -3.02 -13.85 -11.71 -28.93 -23.80 -14.17 -12.35 -34.54

Mulliken Voogd22 Dreiding41 -3.40 -17.68 -14.05 -41.45 -31.18 -17.70 -14.72 -59.43

Figure 3. R-Glycine hydrogen-bonded bilayers, produced using Cerius2.62 Intermolecular interactions, as described in Table 4, are marked.

cant. However, the charge distributions by No and Voogd do not differ considerably; therefore, this cannot be concluded. The model based on the Dreiding force field resulted in a shape, which is very similar to the one based on the Scheraga force field. Because we have already established that the Dreiding model is inappropriate, one may find this result quite surprising. However, the morphological importance of a face is determined by its relative rather than its absolute attachment energy; therefore, when the relative attachment energies are similar, the resulting morphologies will be similar, regardless of the total suitability of the model (Figure 5). Shapes similar to the ones by Boek et al.12 were obtained when the energies were calculated using the Scheraga force field without any charge description. Shapes similar to the ones reported by BerkovitchYellin10a could not be reproduced. When her calculated



Figure 5. Predictions of R-glycine by the attachment energy model. (a) Shape according to attachment energies reported by Berkovitch-Yellin. (b) Shape according to attachment energies reported by Boek et al., model 1. (c) Shape according to the attachment reported by Boek et al., model 3. (d) Shape according to the AE model, using Scheraga force field and charges due to No. (e) Shape according to the AE model, using Dreiding force field and charges due to Voogd. Figure 4. Prediction of R-glycine by the BFDH model in two projections. The external morphology is described along with two unit cells to allow visualization of the packing configuration.

values are used to construct a shape, the resulting shape (shown in Figure 5a) is different than the one reported by her. This was also noted by Boek et al.12 BCF Model Accounting for Solute-Solvent Surface Free Energies. Li and Rodrı´guez-Hornedo68 measured the growth rates of individual faces of R-glycine crystals grown from water and concluded that the growth of glycine from aqueous solution could be explained by the screw dislocation mechanism. Strong intermolecular interactions, “bonds”, were calculated using the Scheraga force field and charges due to No et al. The calculation was performed by means of the program Habit45 and reported in Table 4. The interactions are similar to those calculated by Berkovitch-Yellin10a where the electrostatic term was represented by a multipolar expansion, derived from deformation electron density maps, and to those of Boek12 model 2, which is based on the same electrostatic description. Interaction a is a very strong hydrogen bond interaction between two monomer units in the cell, M(1) and

M(3). The strong nature of this interaction suggests that the growth unit is a dimer, as determined by Gidalevitz et al.55 Therefore, we created a unit cell with two growth units, D(1) comprised of monomers M(1) and M(3) and D(2) comprised of monomers M(2) and M(4). The units are centered at D(1) x, y, z and D(2) x + 1/2, -y + 1/2, z + 1/2. Intermolecular interactions for this new cell were calculated, and they are reported in Table 5. Crystal and kink surface free energies were calculated, and they are reported along with kink surface area and solvent surface free energy in Table 6. High kink energies suggest that crystal faces grow by the screw dislocation mechanism when R-glycine is crystallized from water, as was concluded empirically by Li and Rodrı´guez-Hornedo.68 Application of our shape prediction model for the case of a dimer growth unit reveals that face {020} contains both bond chains, a and b, and each of the following faces: {110}, {120}, {011}, and {021}, contains one bond chain. The predicted shape (Figure 6b) is an elongated prism, which resembles the experimental grown shape (Figure 6a) and includes all of the external faces, except for {120}. A better prediction can be made if one selects the monomer as the growth unit thus considering the

232



Table 4. Intermolecular Interactions for r-Glycine interaction energy (kcal/mol) bond type a b c d

description

attractive

repulsive

Coulombic

total

M(1)-M(3) M(1)-M(1, 001); M(1)-M(1, 001 h) M(1)-M(3, 100) M(1)-M(2); M(1)-M(2, 1 h 01 h)

-6.49 -10.43

4.60 10.59

-9.51 -6.56

-11.401 -6.399

-3.06 -1.98

1.79 1.26

-6.33 3.65

-7.591 2.934

Table 5. Intermolecular Interactions for r-Glycine Based on a Dimer Growth Unit interaction energy (kcal/mol) bond type a b

description

attractive

repulsive

Coulombic

total

M(1)-M(1, 001 h ); M(1)-M(1, 001) M(1)-M(1, 100); M(1)-M(1, 1 h 00)

-22.16

21.67

-10.09

-10.582

-16.05

12.50

-6.56

-10.101

Table 6. Solute, Solvent, and Kink Surface Free Energies for r-Glycine Based on a Dimer Growth Unit bond chain a b

φcryst (kcal/mol)

Akink (Å2/kink)

γcryst (erg/cm2)

γsolv (erg/cm2)

γsolv h (erg/cm2)

γkink (erg/cm2)

φkink (kcal/mol)

10.582 10.101

28.360 30.358

259.5 231.4

72.800 72.800

58.984 58.984

273.3 245.2

11.15 10.71

interactions in Table 4 but excludes interaction a, which has an intramolecular character. This approach is similar to selecting the dimer as the growth unit in the sense that interaction a is considered an intramolecular interaction, but the nature of the unit cell (with four units) and of other interactions stay unchanged. Because interaction d is a repulsion force overall, it cannot form a bond chain and it is irrelevant for the purpose of shape calculation. This calculation was therefore based on merely two strong PBCs: b and c. Faces {111} and {021 h }, which are not parallel to a PBC, are ruled out, and face {101 h }, which is parallel only to PBC a, is not included in the final shape, as this strong interaction is treated as an intramolecular interaction. Faces that contain strong PBCs are {020}, {011}, {110}, and {120}. Calculated solute, solvent, and kink surface free energies are reported in Table 7, and the predicted shape is presented in Figure 6c. When R-glycine is crystallized out of water, it exhibits bipyramidal morphology composed of large {110} and {011} faces and smaller {010} faces.10a,12,69 Boek et al.12 reported in addition to these faces the appearance of {120}, which accounts for a small part of the external side morphology. Kang et al.,69 who also crystallized R-glycine, reported {110} and {020} as side faces for growth out of aqueous solution and {110} and {120} as side faces for growth on a 100% OH self-assembled monolayer (SAM) surface. Our predicted shape, using the second approach, is of an elongated box with tilted edges that resembles a coffin (Figure 6c). External faces include {011}, {020}, {110}, and {120}. There is a strong similarity between the predicted shape and the actual grown shape, as seen in Figure 6a,c. The best calculated shape is obtained using a hybrid monomer/dimer growth unit. It remains an open question how to arrive at the best choice of growth unit in a predictive fashion without the prior knowledge of the observed shape. The best that can be achieved at the

present time in a predictive mode is to select the dimer growth unit. L-Alanine L-Alanine (H2NCH(CH3)COOH), like R-glycine, belongs to the R-amino acids and crystallizes as zwitterions. It is used by the body as a source of energy for muscle tissue, the brain, and the central nervous system.19 Because it is the smallest amino acid with a chiral R-carbon, L-alanine is used to model more complex residues in proteins and peptides.21 L-Alanine crystallizes in the orthorhombic space group P212121, with four molecules per unit cell.71 Unit cell dimensions are a ) 6.025 Å, b ) 12.343 Å, and c ) 5.784 Å. General position coordinates are M(1) x, y, z; M(2) -x + 1/2, -y, z + 1/2; M(3) x + 1/2, -y + 1/2, -z; and M(4) -x, y + 1/2, -z + 1/2, as can be seen in Figure 2. The extinction conditions for the space group are65 h00, h ) 2n; 0k0, k ) 2n; and 00l, l ) 2n. The calculated center of mass is at xCM ) 0.5278, yCM ) 0.1331, and zCM ) 0.4671. To the best of our knowledge, the shape of L-alanine, when crystallized from the vapor, is not reported in the literature. The BFDH rule is used as the first step in the process of face selection. Using the BFDH model to estimate the shape yields the results shown in Figure 7. Attachment Energy Model. Attachment energies were calculated using the program Cerius2 66 with the selected force fields of Scheraga40 and Dreiding41 (reported in Table 8). Partial charges used are due to No et al.,20b,c who carried out ab initio calculation using the 6-31G** basis set and PD method to describe the charge, and due to Voogd et al.,22 who carried out ab initio calculations using the minimal basis set and Mulliken population to describe the charge; both are reported in Table 2. In an effort to test the importance of charge description and force field selection, we calculated the attach-



Figure 8. Morphologies of L-alanine predicted by the attachment energy model. (a) Scheraga force field with charges due to No. (b) Scheraga force field with charges due to Voogd. (c) Dreiding force field with charges due to Voogd. (d) Lifson force field with charges due to No.

Figure 6. Reported and predicted morphologies for R-glycine crystallized out of water. (a) R-Glycine grown out of aqueous solution, from Boek et al.12 (b) Prediction of R-glycine grown out of aqueous solution based on a dimer growth unit. (c) Prediction of R-glycine grown out of aqueous solution based on a modified monomer growth unit.

Figure 7. Prediction of L-alanine by the BFDH model.

ment and lattice energies using the Lifson Hagler Dauber28e force field and additional force fields. The Lifson Hagler Dauber force field was created for carboxylic acids and amides and was not expected to be a good choice for the calculation of the interactions of L-alanine, since amino acids were not included in its development. First, the original set of point charges provided by the force field was assigned to the model. The resulting attachment energies of several faces were positive, which implied that the force field was unsuitable for modeling the system, or that the charge was not properly represented, or both. Next, point charges by No and by Voogd were assigned and used together with the force field. This resulted in overestimated negative lattice and attachment energies, as can be seen in Table 8, which led us to the conclusion that the Lifson force field is unsuitable for modeling this system. The

Scheraga force field is the only one that reproduced the lattice energy in accordance with the sublimation energy; therefore, it is the most appropriate of the three. The shapes based on attachment energy calculations with different force fields are presented in Figure 8, and they are all elongated shapes exhibiting {110} and {020} as the vertical faces and {011} as the top face. Observing the similarity of the shapes, one can mistakenly conclude that any of the selected force fields is a suitable selection for this system. However, as was noted for R-glycine, the attachment energy model is not sensitive to the absolute value of the energies but only to their relative values. Therefore, one cannot decide, merely by observing a shape, that a certain potential is adequate/ inadequate. This lack of sensitivity of the attachment energy model to the applied force filed was also noted by Brunsteiner and Price.67 BCF Model Accounting for Solute-Solvent Surface Free Energies. Intermolecular interactions were calculated using the Scheraga force field and charges due to No et al. There are two strong intermolecular interactions (reported in Table 9). The growth unit initially selected was a monomer, and because the packing orientation and the interactions did not suggest a different selection, it was used as the basic building unit for the model. High kink energies (Table 10) suggest that crystal faces grow by the screw dislocation model when crystallized from water, as was concluded by Lechuga-Ballesteros and Rodrı´guez-Hornedo72b on different grounds. Their study provides a detailed morphology of L-alanine crystallized from water, along with measured growth rates of the appearing faces. These growth rates indicate that for low supersaturation levels, faces {010} and {120} grow at approximately the same growth rate, while face {011} grows at twice this rate. The growth rates describe a prismatic shape, elongated along the c-axis and bounded by {120}, {010}, and {110} faces, as sketched in their article. Faces {011} are reported as the top faces. Our calculations produced similar results; however, we found a difference in the assignment of side faces. To test our calculations, we crystallized L-alanine out of water. Experiments. We grew the crystals in 10 mL glass tubes, immersed in a thermal bath. The materials used were L-alanine, 98% pure, by Sigma-Aldrich and deionized water, treated with a NANOpure infinity system. Crystallization was achieved by slow cooling. The solution was initially heated to dissolve all solids and

234



Table 7. Solute, Solvent, and Kink Surface Free Energies for r-Glycine Based on a Monomer Growth Unit bond chain

φcryst (kcal/mol)

Akink (Å2/kink)

γcryst (erg/cm2)

γsolv (erg/cm2)

γsolv h (erg/cm2)

γkink (erg/cm2)

φkink (kcal/mol)

6.56 6.33

14.180 16.503

321.7 266.7

72.800 72.800

58.984 58.984

335.503 280.529

6.85 6.66

b c

Table 8. Calculated Attachment Energies for Selected Faces of L-Alanine (kcal/mol) face

dhkl

this work PD No17b,c

Mulliken Voogd19

Mulliken Voogd19

PD No17b,c

Mulliken Voogd19

Scheraga40 -10.11 -8.78 -17.52 -11.93 -19.62 -24.75 -24.74 -38.26

Scheraga40 -9.31 -7.99 -18.89 -10.99 -20.73 -26.07 -26.15 -39.50

Dreiding41 -12.07 -9.26 -29.16 -13.35 -31.22 -40.07 -39.81 -56.07

Lifson28e -20.53 -14.87 -41.54 -22.85 -45.84 -56.71 -56.60 -84.56

Lifson28e -19.62 -14.96 -43.06 -22.18 -46.99 -58.55 -58.33 -87.12

charge density force field 020 110 011 120 021 101 111 Elatt

6.172 5.414 5.238 4.311 4.220 4.173 3.953

Table 9. Intermolecular Interactions for L-Alanine interaction energy (kcal/mol) bond type a b

description

attractive

repulsive

Coulombic

total

M(1)-M(3, 1 h 01); M(1)-M(3, 001) M(1)-M(1, 001); M(1)-M(1, 001 h)

-8.45

7.12

-7.75

-9.082

-8.58

8.35

-6.85

-7.095

Table 10: Solute, Solvent, and Kink Surface Free Energies for L-Alanine bond chain a b

(kcal/mol)

Akink (Å2/kink)

γcryst (erg/cm2)

γsolv (erg/cm2)

γsolv h (erg/cm2)

γkink (erg/cm2)

φkink (kcal/mol)

7.75 6.85

25.69 18.63

209.7 255.7

72.8 72.8

58.98 58.98

223.55 269.54

8.26 7.22

φcryst

then slowly cooled, at a very low cooling rate of 0.1-0.2 °C per day. Higher cooling rates produced crystals that were not well-faceted and resulted in the occurrence of spontaneous nucleation. The solution was seeded with a few seeds just prior to entering the metastable zone width. The seeding, combined with temperature control, allowed the growth process to progress while maintaining the solution within the metastable zone. Crystal faces were indexed using a Nonius KappaCCD single crystal X-ray diffractometer equipped with a CCD area detector and a video microscope. The predicted shape of L-alanine grown from aqueous solution, along with the experimentally grown shape, is presented in Figure 9. The prediction compares well with the actual grown shape, except for face {110}, which was not predicted but observed. However, face {110} has low morphological importance, as can be seen in Figure 9d, and its absence does not significantly detract from the generally good prediction. Conclusions We have explored models to predict the shapes of amino acids grown from solution. Amino acids that crystallize in the form of zwitterions are involved in a process of proton transfer. Empirical evidence as well as thermodynamic calculations suggest that for glycine, and possibly for other amino acids, the transfer process occurs in the gas phase. Glycine in solution and in the solid exists in the zwitterionic form. Therefore, despite the large values reported for proton transfer energies,

in the order of sublimation enthalpies, it should have no effect on the growth of amino acid crystals from solution. A reliable description of net atomic charges is important for the calculation of electrostatic interactions of polar molecules in general and of amino acids in particular. Whenever available, ab initio calculated charges by advanced basis sets provide the most accurate description. Otherwise, one can calculate charges by molecular orbital methods, which allow direct calculation based on atomic and orbital characteristics. The selection of a force field should be in accordance with the modeling objective. Although many crystals can be modeled using several force fields, not all force fields provide a good representation of the intermolecular interactions. The fitness of a force field can be assessed by the calculated lattice energy, which should approximately equal the sublimation enthalpy. Determination of the growth unit is not always obvious, and different growth units are often suggested for the same system. The selection of a growth unit influences the crystal shape since it implies that certain interactions are intramolecular, thus excluded from the model, and others intermolecular, thus included in the model. Even though the structure of the growth unit is important for modeling purposes, some studies tend to overlook the importance of its determination and others specify a certain growth unit based on intuition. A graphical observation of the orientation of the molecules in the solid state as well as analysis of the intermolecular interactions were suggested as guiding rules for the



solvent. Additionally, it would be useful to expand the model to account for other external factors such as cosolvents and impurities. Acknowledgment. We are grateful to the sponsors of the Process Design and Control Center at the University of Massachusetts, Amherst. V.B.-L. acknowledges a Fellowship from the LIFE foundation and thanks Prof. L. Leiserowitz for valuable advice. We thank the X-ray Structural Characterization Laboratory at the Chemistry Department supported by the University of Massachusetts and National Science Foundation (Grant CHE-9974648) and Dr. A. Chandrasekaran for data collection for L-alanine. Nomenclature

Figure 9. Shape of L-alanine when grown from aqueous solution. (a) Prediction of L-alanine using our solvent-solute model. (b) The same prediction as in panel a from a different projection. (c) A crystal of L-alanine grown from aqueous solution in our lab. (d) Faces of L-alanine, as indexed by X-ray diffractometery.

determination of a growth unit. Whenever there is doubt, several growth units need to be tested and the selection is not always obvious or unique. For instance, for the R-glycine system, the best results were achieved by using an intermediate approach between monomer and dimer growth units. A model for predicting the crystal shape of polar and nonpolar molecules was presented. Starting with a geometrical crystal model and a list of ∼50 faces dictated by the BFDH rule, it provides rules and an algorithm to reach a final prediction. The methodology is based on a detailed BCF growth mechanism and requires calculations of kink and edge free energies at the solute-solvent interface. It includes special considerations for polar systems. R-Glycine and L-alanine were selected as model compounds for application of the shape prediction model. The predicted shape of R-glycine when grown from aqueous solution resembles a coffin, and it compares well with the experimentally grown shape by Boek et al.12 The predicted shape of the shape of L-alanine compares well with our experimentally grown shape. Future efforts could be usefully directed toward more challenging systems such as amino acids with significant hydrophobic, polar, and charged side chains, crystallized from various organic/aqueous solvents, to test and extend the method to address a wider variety of important forces that may exist between solute and

ahkl, intermolecular distance in the direction of the step (Å) aphkl, the distance the step is propagated by a adding a monolayer to it (Å) C gneu , heat capacity of the gas phase neutral form (energy/mol) C g,x zwi , heat capacity of the gas phase zwitterion whose geometry is the same as that of the zwitterion in the molecular crystal (energy/mol) C gzwi , heat capacity of the gas phase zwitterionic form (energy/mol) C szwi , heat capacity of the solid phase zwitterionic form (energy/mol) dhkl, interplanar spacing (Å) ∆E ab err,neu , the energy difference between the neutral molecule in its equilibrium conformation in the gas at 0 K and the ab initio minimum energy conformation at 0 K (energy/mol) ∆E ab err,zwi , the energy difference between the zwitterion in its equilibrium conformation in the gas phase and the ab initio minimum energy conformation at 0 K (energy/mol) ∆Ept, proton transfer energy at 0 K (energy/mol) ∆Erelax, the relaxation energy when the gas phase zwitterions, with the most stable geometry, are changed into the zwitterions with the same geometry as in a molecular crystal (energy/mol) Elatt, the energy of formation of a crystal from the isolated (gas phase) molecules (energy/mol) hhkl, step height (Å) H g, enthalpy of the gas (energy/mol) H cr, enthalpy of the crystal (energy/mol) ∆Hsub, sublimation enthalpy (energy/mol) P, vapor pressure (Pa) Pref, reference vapor pressure (Pa) qi,qj, partial charge placed on atoms i, j, respectively r, interatomic distance (Å) R, ideal gas constant Rhkl, rate of growth in the direction normal to a surface (Å/s) T, temperature (K) Tref, reference temperature (K) U gvibr , U cr vibr , the additional energy of vibration at temperature T (energy/mol) U gel , U cr el , the total (fixed nuclei) molecular energy (energy/mol) vhkl, lateral step velocity (Å/s)

236


vkij, the interaction energy between atom i in the central molecule and atom j in the kth surrounding molecule (energy/mol) WA, work of adhesion X 0hkl , the average distance between the kinks on a step coinciding with a PBC (Å) yhkl, distance between steps (Å) γl, interfacial free energy (surface tension) of the liquid (energy/mol) γdl , the dispersion force part of the surface free energy of the liquid (energy/mol) γhl , the hydrogen bond part of the surface free energy of the liquid (energy/mol) γs, interfacial free energy (surface tension) of the solid (energy/mol) γls, interfacial free energy at the interface between the liquid and the solid (energy/mol) g0, cr 0 , zero point vibrational energy (energy/mol) φkink hkl , the energy required to form a kink (energy/ mol)


(21) (22) (23) (24)

(25)

(26)

(27) (28)

References (1) Valder, C.; Merrifield, D. Pharmaceutical Technology. SmithKline Beecham R&D News 1996, 32, 1. (2) Shekunov, B. Y.; York, P. J. Cryst. Growth 2000, 211, 122136. (3) Fagan, P. F.; Hammond, R. B.; Roberts, K. J. In Crystal Growth of Organic Materials; Myerson, A. S., Green, D. H., Meenan, P., Eds.; ACS: Washington, DC, 1996; p 22-27. (4) Brittain, H. G. Polymorphism in Pharmaceutical Solids; Marcel Dekker: NY, 1999. (5) Bernstein, J. J. Phys. D: Appl. Phys. 1993, 26, B66-B76. (6) Friedel, M. G. Bull. Soc. Fr. Mineral. 1907, 30, 326-455. (7) Donnay, J. D. H.; Harker, D. Am. Mineral. 1937, 22, 446467. (8) Wells, A. F. Philos. Mag. 1946, 37, 184-199. (b) Wells, A. F. Philos. Mag. 1946, 37, 217-236. (9) Hartman, P.; Perdok, W. G. Acta Crystallogr. 1955, 8, 4952. (b) Hartman, P. In Physics and Chemistry of the Organic Solid State; Fox, D., Labes, M. M., Weissberger, A., Eds.; Wiley: New York, 1963; p 369-409. (c) Hartman, P. In Crystal Growth: An Introduction; Harman, P., Ed.; NorthHolland: Amsterdam, 1973; p 367. (d) Hartman, P.; Bennema, P. J. Cryst. Growth 1980, 49, 145-156. (10) Berkovitch-Yellin, Z. J. Am. Chem. Soc. 1985, 107, 82398253. (b) Berkovitch-Yellin, Z.; van Mil, J.; Addadi, L.; Idelson, M.; Lahav, M.; Leiserowitz, L. J. Am. Chem. Soc. 1985, 107, 3111-3122. (11) Lahav, M.; Leiserowitz, L. Chem. Eng. Sci. 2001, 56, 22452253. (12) Boek, E. S.; Feil, D.; Briels, W. L.; Bennema, P. J. Cryst. Growth 1991, 114, 389-410. (13) Docherty, R.; Roberts, K. J. J. Cryst. Growth 1988, 88, 159168. (14) Winn, D.; Doherty, M. F. AIChE J. 1998, 44, 2501-2514. (15) Liu, X.-Y.; Bennema, P. Phys. Rev. E 1993, 48, 2006-2015. (b) Liu, X.-Y.; Bennema, P. J. Chem. Phys. 1993, 98, 58635872. (c) Liu, X.-Y.; Bennema, P. Phys. Rev. B 1994, 49, 765-775. (d) Liu, X.-Y.; Boek, E. S.; Briels, W. J.; Bennema, P. Nature 1995, 374, 342-345. (e) Liu, X.-Y.; Bennema, P. Phys. Rev. B 1996, 53, 2314-2321. (f) Liu, X.-Y.; Bennema, P. J. Cryst. Growth 1996, 166, 112-123. (16) Bisker-Leib, V.; Doherty, M. F. Cryst. Growth Des. 2001, 1 (6), 455-461. (17) Turner, A. J.; Whittle, S. R. Biochem. J. 1983, 209, 29-41. (18) Weber, H.-P.; Craven, B. M.; McMullan, R. K. Acta Crystallogr. 1983, B39, 360-366. (19) Wolff, M., Ed. In Burger’s Medicinal Chemistry, 4th ed; Wiley: New York, 1979. (20) No, K. T.; Grant, J. A.; Scheraga, H. A. J. Phys. Chem. 1990, 94, 4732-4739. (b) No, K. T.; Grant, J. A.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1990, 94, 4740-4746. (c)

(29)

(30)

(31)

(32) (33) (34) (35) (36) (37)

(38)

(39)

(40) (41) (42) (43)

No, K. T.; Cho, K. H.; Kwon, O. Y.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1994, 98, 10742-10749. (d) Kwon, O. Y.; Kim, S. Y.; No, K. T.; Kang, Y. K.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1996, 100, 17670-17677. Csa´sza´r, A. G.; Perczel, A. Prog. Biophys. Mol. Biol. 1999, 71, 243-309. Voogd, J.; Derissen, J. L.; van Duijneveldt, F. B. J. Am. Chem. Soc. 1981, 103, 7701-7706. Brown, R. D.; Godfrey, P. D.; Storey, J. W. V.; Bassez, M.P. J. Chem. Soc. Chem. Commun. 1978, 547-548. Bernstein, J.; Etter, M. C.; Leiserowitz, L. In Structure Correlation; Bu¨rgi, H.-B., Dunitz, J. D., Ed.; VCH: Weinheim, 1994; Vol. 2, pp 431-507. Addadi, L.; Berkovitch-Yellin, Z.; Weissbuch, I.; van Mil, J.; Shimon, L. J. W.; Lahav, M.; Leiserowitz, L. Angew. Chem., Int. Ed. Engl. 1985, 24, 466-485. Weissbuch, I.; Popovitz-Biro, R.; Leiserowitz, L.; Lahav, M. In The Lock-and-Key Principle; Behr, J.-P., Ed.; John Wiley & Sons Ltd.: New York, 1994; pp 173-246. Gavish, M.; Wang, J.-L.; Eisenstein, M.; Lahav, M.; Leiserowitz, L. Science 1992, 256, 815-818. Hagler, A. T.; Huler, E.; Lifson, S. J. Am. Chem. Soc. 1974, 97, 5319-5326. (b) Hagler, A. T.; Lifson, S. J. Am. Chem. Soc. 1974, 97, 5327-5335. (c) Hagler, A. T.; Leiserowitz, L.; Tuval, M. J. Am. Chem. Soc. 1976, 98, 4600. (d) Bernstein, J.; Hagler, A. T. J. Am. Chem. Soc. 1978, 100, 673. (e) Lifson, S.; Hagler, A. T.; Dauber, P. J. Am. Chem. Soc. 1979, 101, 5111-5121. (f) Lifson, S.; Hagler, A. T.; Dauber, P. J. Am. Chem. Soc. 1979, 101, 5122-5130. (g) Hagler, A. T.; Dauber, P.; Lifson, S. J. Am. Chem. Soc. 1979, 101, 5131-5141. Gilli, G. In Fundamentals of Crystallography; Giacovazzo, C., Ed.; Oxford University Press: New York, 1995; pp 465534. Docherty, R. In Crystal Growth of Organic Materials; Myerson, A. S., Green, D. H., Meenan, P., Eds.; ACS: Washington, DC, 1996; p 2-14. Chickos J. S. Heats of Sublimation. In Molecular Structure and Energetics; Liebman, J. F., Greenberg, A., Eds.; Vol. 2, Physical Measurements; New York, 1987; pp 67-150. Mirsky, K. Acta Crystallogr. 1976, A32, 199-207. Takagi, S.; Chihara, H.; Seki, S. Bull. Chem. Soc. Jpn. 1959, 32, 84-88. Born, M. Z. Phys. 1920, 1, 45. Rae, A. I. M.; Mason, R. Proc. R. Soc. London, Ser. A 1968, A304, 487-499. Tortonda, F. R.; Pascual-Almir, J. L.; Silla, E.; Tuno´n, I. Chem. Phys. Lett. 1996, 260, 21-26. Kitaigorodskii, A. I. Organic Crystal Chemistry; Consultants Bureau: New York, 1961. (b) Kitaigorodskii, A. I. Acta Crystallogr. 1965, 18, 585. (c) Kitaigorodskii, A. I.; Mirskaya, K. V.; Tovbis, A. B. Sov. Phys. Cryst. 1968, 13, 176-180. (d) Kitaigorodskii, A. I. Advances in Structure Research by Diffraction Methods; Brunswick: 1970; Vol. 3. (e) Perstin, A. J.; Kitaigorodskii, A. I. The Atom-Atom Potential Methodol. Application to Organic Molecular Solids; SpringerVerlag: New York, 1985; p 79. Williams, D. E. J. Chem. Phys. 1965, 43, 4424. (b) Williams, D. E. Science 1965, 147, 605. (c) Williams, D. E. J. Chem. Phys. 1966, 45 (10), 3770. (d) William, D. E.; Cox, S. R. Acta Crystallogr., Sect. B 1984, 40, 404. (e) Williams, D. E.; Houpt, D. J. Acta Crystallogr. Sect. B 1986, 42, 286. Yan, J. F.; Momany, F. A.; Hoffman, R.; Scheraga, H. A. J. Phys. Chem. 1970, 74, 420-433. (b) Momany, F. A.; McGuire, R. F.; Yan, J. F.; Scheraga, H. A. J. Phys. Chem. 1970, 74, 2424-2438. (c) McGuire, R. F.; Momany, F. A.; Scheraga, H. A. J. Phys. Chem. 1972, 76, 375. (d) Momany, F. A.; Carruthers, L. M.; McGuire, R. F.; Scheraga, H. A. J. Phys. Chem. 1974, 78 (16), 1595-1620. (e) Momany, F. A.; Carruthers, L. M.; Scheraga, H. A. J. Phys. Chem. 1974, 78 (16), 1621-1630. Nemethy, G.; Pottle, M. S.; Scheraga, H. A. J. Phys. Chem. 1983, 87, 1883-1887. Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III. J. Phys. Chem. 1990, 94, 8897-8909. Chipot, C.; Angyan, J. G.; Maigret, B.; Scheraga, H. A. J. Phys. Chem. 1993, 97, 9788-9796. Scott, R. A.; Scheraga, H. A. J. Am. Phys. 1966, 45, 2091.

Modeling Crystal Shape of Polar Organic Materials (44) Chin, D. N.; Zerkowski, J. A.; Macdonald, J. C.; Whitesides, G. M. In Organized Molecular Assemblies in the Solid State; Whitesell, J. K., Ed.; John Wiley & Sons: New York, 1999; pp 185-253. (45) Clydesdale, G.; Docherty, R.; Roberts, K. J. Comput. Phys. Commun. 1991, 64, 311-326. (46) Berkovitch-Yellin, Z.; Leiserowitz, L. J. Am. Chem. Soc. 1982, 104, 4052-4064. (47) Gasteiger J.; Marsili, M. Tetrahedron 1980, 36, 3219-3288. (48) Gavezzotti, A. In Theoretical Aspects and Computer Modeling of the Molecular Solid State; Gavezzotti, A., Ed.; John Wiley and Sons: New York, 1997; pp 1-30. (b) Gavezzotti, A. In Structure Correlation; Burgi, H.-B., Dunitz, J. D., Ed.; VCH: Weinheim, 1994; Vol. 2, pp 510-542. (c) Gavezzotti, A.; Filippini, G.; Kroon, J.; van Eijck, B. P.; Klewinghaus, P. Chem. Eur. J. 1997, 3 (6), 893-899. (49) Hinze, J.; Jaffe, H. H. J. Phys. Chem. 1963, 67, 1501-1505. (50) Svec, H. J.; Clyde, D. D. J. Chem. Eng. Data 1965, 10 (2), 151-152. (51) Davey, R. J.; Blagden, N.; Rignini, S.; Alison, H.; Quayle, M. J.; Fuller, S. Cryst. Growth Des. 2001, 1 (1), 59-65. (52) Grimbergen, R. F. P.; Bennema, P. In Crystal Growth of Organic Materials; ACS Conference Proceedings Series; American Chemical Society: Washington, DC, 1996; pp 2835. (53) Chang, Y. C.; Mayerson, A. S. AIChE 1986, 32 (9), 15671569. (54) Myerson, A. S.; Lo, P. Y. J. Cryst. Growth 1990, 99, 10481052. (55) Gidalevitz, D.; Feidenhans’l, R.; Matlis, S.; Smilgies, D.-M.; Christensen, M. J.; Leiserowitz, L. Angew. Chem., Int. Ed. Engl. 1997, 36, 955-959. (56) Sarma, J. A. R. P.; Desiraju, G. R. Cryst. Growth Des. 2002, 2 (2), 93-100. (57) Burton, W. K.; Cabrera, N.; Frank, F. C. Philos. Trans. R. Soc. 1951, A243, 299.

Crystal Growth & Design, Vol. 3, No. 2, 2003 237 (58) Chernov, A. A. Modern Crystallography III; SpringerVerlag: New York, 1984. (b) Chernov, A. A.; Nishinaga, T. In Morphology of Crystals, Part A; Sunagawa, I., Ed.; Terra Scientific: Tokyo, 1987. (59) Land, T. A.; Malkin, A. J.; Kutznesov, Y. G.; McPherson, A.; De Yoreo, J. J. J. Cryst. Growth 1996, 166, 893-899. (60) Geil, P. Polymer Single Crystals; Interscience: 1963. (61) Kaichev, R. In Growth of Crystals 3; Shubnikov, A. V., Sheftal, N. N., Eds.; Consultants Bureau: New York, 1962. (62) Girifalco, L. A.; Good, R. J. J. Phys. Chem. 1957, 61, 904909. (63) Fowkes, F. M. In Contact Angle, Wettability and Adhesion; Advances in Chemistry 43; Gould, R. F., Ed.; American Chemical Society: Washington, DC, 1964; Chapter 6, pp 99111. (64) Jonsson, P. G.; Kvick, A. Acta Crystallogr. 1972, B28, 18271833. (65) International Tables of Crystallography; Reodel: Dordrecht, 1983; Vol. A. (66) Cerius2 Version 4.2 MatSci; Molecular Simulation Inc.: San Diego, CA, 2000. (67) Burnsteiner, M.; Price, S. L. Cryst. Growth Des. 2001, 1 (6), 447-453. (68) Li, L.; Rodrı´guez-Hornedo, N. J. Cryst. Growth 1992, 121, 33-38. (69) Kang, J. F.; Zaccaro, J.; Ulman, A.; Myerson, A. Langmuir 2000, 16, 3791-3796. (70) Lin, C. H.; Gabas, N.; Canselier, J. P.; Pe`pe, G. J. Cryst. Growth 1998, 191, 791-802. (71) Lehman, M. S.; Koetzle, T. F.; Hamilton, W. C. J. Am. Chem. Soc. 1972, 94, 2657-2660. (72) Lechuga-Ballesteros, D.; Rodrı´guez-Hornedo, N. J. Colloid Interface Sci. 1993, 157, 147-153. (b) Lechuga-Ballesteros, D.; Rodrı´guez-Hornedo, N. Int. J. Pharm. 1995, 115, 139149. (73) Cheetham, A. K. UCSB, private communication, 2001.

CG025538Q

Modeling Crystal Shape of Polar Organic Materials - American

Recommend Documents