Elastic Properties, Thermal Expansion, and Polymorphism of

May 27, 2010 - The thermal expansion is anisotropic, with a minimum of the longitudinal .... coefficients of the tensor of linear thermal expansion (R...
0 downloads 0 Views 2MB Size
DOI: 10.1021/cg100241c

Elastic Properties, Thermal Expansion, and Polymorphism of Acetylsalicylic Acid

2010, Vol. 10 3132–3140

Johannes D. Bauer, Eiken Hauss€ uhl, Bj€ orn Winkler,* and Dirk Arbeck Goethe Universit€ at Frankfurt, Institut f€ ur Geowissenschaften, Abt. Kristallographie, Altenh€ oferallee 1, 60438 Frankfurt am Main, Germany

Victor Milman and Struan Robertson Accelrys Inc., 334 Cambridge Science Park, Cambridge CB4 0WN, United Kingdom Received February 22, 2010; Revised Manuscript Received May 17, 2010

ABSTRACT: High quality single crystals of acetylsalicylic acid (aspirin) with dimensions up to about 30 cm3 were grown from saturated solutions. The complete elastic stiffness tensor was determined by two different ultrasonic techniques. The thermal expansion coefficients were measured by dilatometry. The thermal expansion is anisotropic, with a minimum of the longitudinal effect parallel to the direction [010]. The expansion perpendicular to this direction is larger by a factor of 2. The anisotropy of the longitudinal elastic stiffness is less pronounced; the difference between the minimum and maximum values is approximately 30%. The maximum elastic stiffness is observed along the directions of the hydrogen bonds linking the centrosymmetric acetylsalicylic acid dimers. The experimentally determined elastic stiffness coefficients and the observed morphology have been compared to results from atomistic model calculations. This comparison allows us to assess the reliability of predictions made with respect to a potential polymorph.

Introduction 1

Since its discovery by Hoffmann at the end of the 19th century, acetylsalicylic acid, C9H8O4, has become a well-known active pharmaceutical ingredient, API. Similar to discussions for other APIs, the question of polymorphism, with its associated medical and legal implications, has been investigated in the past few years.2-5 It turned out that acetylsalicylic acid is a nontrivial case, where, in addition to the well-known form I, some crystallization conditions lead to the formation of crystals with stacking faults. Diffraction data from such crystals can be interpreted as indicating the intergrowth of two types of domains,4 the first of which has the well-known structure of form I,6 while the domains of form II have a local structure corresponding to that found by Vishweshwar et al.,2 where molecules are connected by catemeric hydrogen bonds. The discussion on the polymorphism in acetylsalicylic acid was triggered partially by results of static atomistic model calculations, which implied that a polymorph similar to form II is slightly more stable (in the athermal limit) than form I but may be unstable with respect to shear.2,7 The suggested low shear coefficients would provide an explanation for the failure to obtain single phase crystals of form II. However, the reliability with which elastic properties of acetylsalicylic acid can be predicted has not been properly evaluated, as not all elastic stiffness coefficients of acetylsalicylic acid have been determined yet. A part of the elastic coefficients of acetylsalicylic acid have been investigated by Kim et al.8 and Ko et al.9 with Brillouin spectroscopy, and there are significant inconsistencies (up to 15%) between the two partial data sets available. For acetylsalicylic acid form I, Kim et al.8 optimized potential parameters of a model based on pair potentials and *To whom correspondence should be addressed. E-mail: b.winkler@ kristall.uni-frankfurt.de. Telephone: þ49 (0)69 79840107. Fax: þ49 (0)69 79840109. pubs.acs.org/crystal

Published on Web 05/27/2010

achieved a satisfactory description of structural parameters. They compared the computed set of elastic stiffness coefficients to data available at that time, and they concluded that with the exception of c22 the agreement for the main diagonal coefficients was satisfactory. Ouvrard and Price7 found that with their force field a second polymorph had a slightly lower total energy (0.2 kJ/mol) in the athermal limit, but concluded that the compound was “unlikely to be an observable polymorph” due to low shear coefficients. From structural considerations it is not obvious why there should be low shear coefficients in form II, as the number of hydrogen bonds is the same in both structures, and their strengths should be comparable, given that the packing density is rather similar. Relative stabilities of form I and form II have also been investigated with quantum mechanical models.5 The DFT-B3LYP calculations in that study required the use of fixed lattice parameters. This limits the accuracy of these calculations, as the remnant stress will be different in the two polymorphs after the “internal” geometry optimizations. The differences in the lattice energies between form I and form II were ≈2 kJ/mol; the difference between calculations for the same polymorph, with slightly different basis sets, was about 1 kJ/mol. These energy differences are too small to be significant. Elastic stiffness coefficients have not been computed with quantum mechanical methods up to now. These results raise the question on how well the atomistic models actually predict structure-property relations beyond a comparison of structural parameters and total energies in the athermal limit. The next level of sophistication is to compare derivatives of the total energy. Several alternatives exist. Vibrational spectroscopy (mainly Raman and infrared, but also neutron molecular spectroscopy) is a widely used tool for the characterization of organic crystals by determining molecular and lattice vibrations. The experimentally determined frequencies can be compared to second derivatives of the total r 2010 American Chemical Society

Crystal Growth & Design, Vol. 10, No. 7, 2010

Table 1. Parameters and Results of the Crystal Structure Determination of Acetylsalicylic Acid at 100 K and of Force Field Calculationsa exp X-ray 100 K form I a (A˚) b (A˚) c (A˚) β (deg) V (A˚3) F (g/cm3) λ (A˚) monochromator θ-interval (deg) absorption coefficient μ (1/mm) absorption correction no. of reflections no. of independent reflections Rσ (%) Rint (%) R1 (%) wR2 (%)

11.2545(3) 6.5467(1) 11.2581(3) 95.910(2) 825.09(3) 1.450 0.71073 graphite 4.79-29.34 0.116

force field model

exp X-ray

form I

form II

100 K form II2

11.5693 6.2963 11.2924 96.632 817.069 1.4646

12.1395 6.3837 11.2571 110.726 815.905 1.4667

12.095(7) 6.491(4) 11.323(6) 111.509(9) 827.1(8) 1.447

empirical 13150 2077

3133

Crystal Structure In order to ascertain that the products obtained by our crystal growth technique were pure form I crystals, we performed an X-ray single crystal structure determination with a κ -geometry diffractometer (XCalibur, Oxford Diffraction) at a temperature of 100 K. The most important results and experimental parameters are summarized in Table 1. For the structure solution and refinement, the shelx97-package11 was used. The lattice parameters are in good agreement with those obtained by Aubrey-Medendorp et al.12 at 90 K with X-ray diffraction (max. deviations about 0.1%) and those published by Wilson,13 taken from neutron diffraction data collected at 100 K (max. deviations about 0.2%); see Table 2. Figure 1 shows the obtained crystal structure, which corresponds to the well-known form I. The R-values obtained here are in very good agreement with those of Bond et al.,4 which is indicative of the absence of any domains of form II. Determination of Physical Properties

εij ¼ Rij ΔT þ βij ðΔTÞ2 þ γij ðΔTÞ3 þ ...

)

The physical properties obtained here are described in a Cartesian reference system (crystal physical reference system) eBi, which is defined relative to the axes aBi of the conventional monoclinic crystallographic coordinate system and its reciprocal lattice, aBi *, by the vectors eB3 cB, eB2 bB* bB and eB1 = eB2  eB3. Determination of the Thermal Expansion Tensor rij. The coefficients of the tensor of linear thermal expansion (Rij) were determined by measuring the deformation (ε) of a sample (parallel plate) while applying a temperature change ΔT. The terms βij and γij in eq 1 are the higher-order thermal expansion coefficients.14 )

energy with respect to atomic displacements. However, such measurements are mainly local probes of bonds in the long wavelength limit. While the wave vector dependence of high frequency phonons can be measured in principle, these experiments are challenging and tedious.10 In contrast, elastic stiffness coefficients sensitively depend on the spatial arrangement of interatomic interactions, especially in compressible low symmetry compounds. Experimentally determined elastic stiffness coefficients can be compared to the second derivative of the total energy with respect to strain and therefore offer an opportunity for a more stringent evaluation of atomistic models. The drawback of this method is that perfect single crystals are required and measurements can become complex for low symmetry compounds. However, as large crystals of acetylsalicylic acid can be grown from solution, we initiated a study to provide the missing data and evaluate the reliability of previous results by determining the elastic behavior of acetylsalicylic acid with two different ultrasound techniques (one using plane waves and parallel plates, while the second method was resonant ultrasound spectroscopy) at room temperature and ambient pressure. In addition, the thermal expansion was studied by dilatometry. The results obtained here are discussed with respect to the reliability of density functional theory-based and force field-based calculations, and the prediction of a shear instability of form II is assessed.

)

Article

ð1Þ

Experimentally determined deformation curves are fitted with a third-order polynomial to yield the coefficients Rij, βij, and γij. Determination of Elastic Stiffness Tensor cij Using the Plane Parallel Plate Ultrasonic Technique. The determination of the elastic stiffness tensor is based on the basic elastodynamic

1780 748

1.37 2.08 3.80 10.99

The structural data of Vishweshwar et al.2 for the so-called form II are also given (see Discussion). All structures are monoclinic with space group symmetry P21/c (IT No. 14) and Z = 4 formula units per unit cell. a

Figure 1. Crystal structure of acetylsalicylic acid.

Table 2. Lattice Parameters of Acetylsalicylic Acid Form I ref

a (A˚)

b (A˚)

c (A˚)

β (deg)

this study, 100 K Aubrey-Medendorp et al.,12 90 K Wilson,13 100 K, neutron-diffr Wheatley,6 ca. 300 K Kim et al.,8 ca. 300 K Wilson,13 300 K, neutron-diffr

11.2545(3) 11.242(7) 11.233(3) 11.446(13) 11.430(1) 11.416(5)

6.5467(1) 6.539(4) 6.540(1) 6.596(6) 6.591(1) 6.598(2)

11.2581(3) 11.245(9) 11.271(3) 11.388(9) 11.395(2) 11.483(5)

95.910(2) 95.9(3) 95.89(2) 95.55(3) 95.68(1) 95.60(3)

3134

Crystal Growth & Design, Vol. 10, No. 7, 2010

Bauer et al.

equations. These relate the propagation velocity of an ultrasonic wave v = fλ to its displacement vector ξk. ð - Fv2 δik þ cijkl gj gl Þξk ¼ 0

ð2Þ

In this equation gj and gl are the components of the propagation directions, δij is the Kronecker symbol, and F is the density. The equation can be solved for the eigenfrequencies for any ξkvector. With the knowledge of the displacement vector, the relationship between the velocity and the cij can be obtained.14 For example, for a longitudinal ultrasonic wave that travels perpendicular to the face of a (010) plate in the case of a monoclinic crystal, one obtains the following: Fv2 ¼ c22

ð3Þ

If a transversal ultrasonic wave with displacement vector ξ 3 eB2 = 0 and propagation direction parallel to eB2 is used for excitation on a sample with the same orientation, the solution of the equation is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 ðc44 þ c66 Þ2 þ 4c246 ð4Þ Fv2 ¼ ðc44 þ c66 Þ ( 2 2 The resonance frequencies can be measured through the excitation of mechanical oscillations in a parallel plate using an ultrasonic transducer. The ultrasonic wave generated at the face of the plate which is in contact with the transducer travels through the crystal and is reflected at the upper face of the plate. It returns a part of its mechanical energy to the transducer via the piezoelectric effect. Resonance frequencies can be detected by measuring the phase angle between the current and the voltage.15 The wave velocity can be obtained from the known thickness of the plate (d) and the difference between two neighboring resonance frequencies Δf (eq 5). The average value of Δf can be obtained from the slope of the linear relationship between the number of resonance modes and the frequency. v ¼ 2dΔf

ð5Þ

Alternatively, one can calculate the wave velocity with two randomly selected resonance modes numbered n and m using the following equation. fn - fm ð6Þ v ¼ 2d n-m Resonant Ultrasound Spectroscopy. Resonant ultrasound spectroscopy (RUS) is based on the determination of the eigenfrequencies of a freely vibrating crystal.16 The spectra are typically obtained in the region between 50 and 2000 kHz. The experimentally determined spectra cannot be directly transformed into the tensor components cij. Instead, the components cij are obtained from the measured resonance frequencies via a nonlinear least-squares refinement procedure in which the observed values are compared to the theoretical ones.17 The theoretical resonance frequencies are obtained by solving a general eigenvalue problem, whose rank is equal to the number of basis functions needed for the description of the components of the displacement vector.16 2 Γa B ¼ ðFω ÞE a B

ð7Þ

The eigenvalues are λ = Fω2, and the eigenvectors, aB, are composed of the expansion coefficients. E is the matrix of the potential energy, and Γ is the matrix of the corresponding

kinetic energy. In the least-squares refinement procedure, the quantity χ (eq 8) is calculated for n eigenmodes with the resonance frequencies fi = ωi/2π and minimized by varying the elastic stiffness coefficients. n X χ ¼ wi ðω2i, cal - ω2i, obs Þ2 ð8Þ i¼1

The wi are the weights calculated by assuming an experimental error of (0.05 kHz for each of the observed resonance frequencies.18 In order to minimize errors due to truncation effects, up to several thousand normalized Legendre polynomials are used in the development of the displacement vector. Computational Details Force field based calculations of the elastic constants were performed using the Forcite program.19 The choice of force field is clearly the principal factor in determining the quality of the predictions, and for molecular crystal systems the terms relating to nonbond interactions are critical to property prediction. Thus, force fields that have been validated against experiment for properties that are critically dependent on nonbond interactions are more likely to generate consistent predictions for other nonbond dependent properties. With this in mind, the force field calculations were performed using the COMPASS20-23 force field, which was developed, primarily, for investigating polymer materials. In brief, valence term parameters were obtained from ab inito calculations adjusted by scaling factors as to give agreement with spectroscopic data. Atomic charges were calculated by the bond increment method, which is based on electrostatic potential calculations and reflects the relative electronegativities between atoms, and van der Waals potential parameters were determined by extensive molecular-dynamics simulation calculations, in which the parameters are adjusted to reproduce condensed phase properties such as mass and cohesive energy densities. In all calculations, Coulomb and van der Waals terms were evaluated using Ewald sums to a precision of ≈ 4 10 -5 kJ/mol. Elastic constants were computed by applying a series of finite strain patterns (consistent with symmetry) to an optimized crystal structure, relaxing the internal degrees of freedom at each strained configuration and obtaining the elastic constant from a linear fit of the calculated stress against strain for the series of strains. In each case, the calculation of elastic constants is preceded by a geometry optimization. During preoptimization, cell parameters can be optimized and separate calculations were performed for each structure in which cell parameters were fixed and variable. When cell parameters were allowed to vary, significant changes were observed, particularly in the case of the acetylsalicylic acid form II structure, and this observation motivated elastic constant calculations for fixed cell geometries. In all cases, the convergence criteria for geometry optimization were a maximum energy change of 8.4  10-5 kJ/mol and a maximum force of 4.2  10-3 kJ/(mol A˚). For each elastic constant, four finite strains where applied with a maximum amplitude of 0.003. In addition, density functional theory based calculations were performed using the CASTEP program.19,24 CASTEP uses pseudopotentials to represent ion-electron interactions; we used on-the-fly generated ultrasoft potentials in this study. The pseudopotential approach makes it possible to represent electronic wave functions in a computationally convenient

Article

Crystal Growth & Design, Vol. 10, No. 7, 2010

plane wave basis set. The energy cutoff of the basis set, one of the most important parameters determining convergence of calculations, was set at 610 eV. Brillouin zone summation was carried out using a 2  4  2 grid of Monkhorst-Pack special points which corresponds to the reciprocal space distances between points of about 0.04 1/A˚. The Wu-Cohen gradient corrected exchange correlation functional25 was employed in all calculations. Our investigations of nearly 100 inorganic compounds showed that this functional produces elastic constants in the best agreement with experiment when compared to standard PBE26 or LDA functionals; it is marginally more accurate than the recent PBEsol functional.27 Calculations were carried out with fixed cell parameters, since it is well-known that the neglect of van der Waals interactions in the standard DFT functionals makes accurate determinations of cell geometries problematic for molecular structures. Geometry optimization of internal degrees of freedom was performed with the convergence criteria of 0.01 eV/A˚ for the highest force component and 0.005 meV for the energy change. Elastic constants were calculated by applying a finite strain to the unit cell and calculating the resultant stress tensor after relaxing internal degrees of freedom in the strained structures.28 The same stringent convergence criteria as for the initial structure optimization were applied. A linear fit of the obtained strain-stress relationships provides second order elastic constants and their statistical uncertainties. Statistical errors were found to be between 0.3 and 1 GPa; hence, the small off-diagonal elements of the calculated cij tensor might not be sufficiently accurate for the purpose of comparison to experimental results.

3135

Figure 2. Crystal of acetylsalicylic acid grown from ethanolic solution by lowering of the temperature. Its thickness is about 5 mm. The bar denotes 10 mm.

Experimental Section Crystal Growth and Sample Preparation. Optical quality crystals were grown from saturated solutions of commercially available acetylsalicylic acid (99%, Alfa Aesar GmbH) in ethanol (pure ethanol, g99.5%, Ph. Eur., Carl Roth GmbH, and denatured ethanol, g99.8% with 1% butanone, Carl Roth GmbH). The growth experiments were carried out in glass vessels in a custom built air thermostat. The temperature was lowered within the range between about 300 and 290 K over periods of 50-60 days. The crystals had sizes between 10  20  25 mm3 and 10  45  75 mm3 (Figure 2). The morphology of the crystals can be described to be between two “extreme” forms, which are both depicted in Figure 3. The main crystal faces of all crystals are {100}, {001}, and {110}, which is in good agreement with the results of Aubrey-Medendorp et al.12 The smaller faces shown in Figure 3 (bottom) are not always present. The measurements of the thermal expansion and the elastic stiffness coefficients were conducted on oriented plane parallel plates. The plates had dimensions of approximately 12  12  5 mm3 with parallelism better than (5 μm ((0.02°). The cuboidshaped sample for the RUS had dimensions of 3.665(4) 4.734(4) 5.063(4) mm3 with its edges parallel to the faces of the crystal physical reference system. The samples were cut by a low speed tungstenwire saw with the aid of a slurry of SiC powder and water. They were further polished using a slurry of corundum powder (grain size 9 μm) and paraffin oil. The deviations from the ideal crystallographic orientation were controlled using an X-ray diffractometer with BraggBrentano geometry (Siemens D500). The deviations were smaller than (0.5°. Density Measurement. The density of acetylsalicylic acid was measured using the buoyancy method in paraffin oil at a temperature of 295.0(5) K. A Sartorius ME235S analytical balance was used. We obtained a density of 1.398(2) g/cm3. This density is in good agreement with the values obtained from the literature. Wheatley6 gives a measured density of 1.40 g/cm3 and a calculated density of 1.398 g/cm3; the calculated density given by

Figure 3. “Extreme” forms of the morphology of acetylsalicylic acid crystals grown from ethanol. Kim et al.8 is 1.401 g/cm3. The calculated X-ray density is often up to 1% larger than the value measured by the buoyancy method, due to defects. The good agreement seen here emphasizes the excellent quality of the crystals grown in the present study. Measurement of Thermal Expansion. The thermal expansion was measured with a commercial push rod dilatometer (Netzsch DIL 402C) equipped with an inductive displacement transducer. Cooling was achieved with a commercial liquid N2 cryostat. During the measurement, the sample environment was flushed with dried helium gas. The coefficients Rii are obtained from samples with surface normals oriented along the principal directions of the crystal physical reference system (eBi), and the coefficient R13 is obtained from measurements on samples with surface normals oriented along (eB1 ( B e 3). Thermal expansion was measured while cooling the samples from 323 to 93 K and reheating them back to room temperature with a temperature increase of 0.5 K/min. In Figure 4 plots of temperature versus strain for the three principal directions of the crystal physical reference system obtained while cooling the samples are depicted. The results were obtained as described above and are given in Table 3; the errors have been obtained by averaging the different measurement runs. Measurement of Elastic Coefficients by Plane Waves Ultrasound Spectroscopy. The coefficients of the elastic stiffness were determined from samples with their surface normals oriented along the

3136

Crystal Growth & Design, Vol. 10, No. 7, 2010

Bauer et al.

Figure 4. Thermal expansion in the main directions of the crystal physical reference system. Table 3. Coefficients of Thermal Expansion of Acetylsalicylic Acid at 293 K

ij

Rij (1  10-6/K)

βij (1  10-9/K)

γij (1  10-12/K)

11 22 33 13

92.6(7) 36.8(1) 80.4(2) 7.7(1)

73(5) -3(1) 78.9(5) 22.4(2)

98(9) -98(4) -91.9(3) -67(4)

Figure 6. Ultrasound resonance spectrum of a plane parallel plate of acetylsalicylic acid with orientation (010) (thickness 5.299(4) mm). For the measurement a longitudinal quartz transducer was used, and its eigenfrequency can bee seen at about 15 MHz.

Figure 5. A sample of acetylsalicylic acid contacted to the ultrasound transducer. directions B e i, (eB1 ( B e 3), and parallel to the natural faces (001), (110), (011), and (112). Three measurements were performed on each of the samples, one with a longitudinal ultrasound transducer (metalized x-cut quartz) and two with a transversal ultrasound transducer (metalized y-cut quartz). The samples were held in contact with the transducer with paraffin oil or resin, the transducer was connected to a network analyzer (Agilent HP 4395A), and the sample holder is depicted in Figure 5. The network analyzer was used to generate an oscillating electric signal, that causes the transducer to perform mechanical oscillations. The analyzer recorded the response signal from the coupled system of the transducer and the sample. The frequency was varied from 4 to 20 MHz in steps of 1 kHz. The temperature during the different measurements was 295.5(10) K. The ultrasonic wave velocity inside the crystal is calculated as described in eq 6 using the lowest- and highest-frequency resonance modes from the impedance spectra. An example spectrum is depicted in Figure 6. Figure 7 depicts the evaluation of the value of Δf from the resonance modes taken from the spectrum of the (010) plate measured with a longitudinal quartz transducer.

Figure 7. Number of resonance modes of the plane parallel plate with orientation (010) plotted against frequency. The coefficients cij of the elastic stiffness tensor were generated from the 27 measured wave velocities using a least-squares algorithm29 and are listed in Table 4. Resonant Ultrasound Spectroscopy. The determination of the elastic properties was also carried out by RUS, using a custom built device (see Figure 8). The sample was gently fixed between two ultrasound transducers. One of the transducers acts as an ultrasound generator, and the other one is employed as an ultrasound detector. Although the force acting on the opposed corners of the cuboid was less than 0.05 N, the acetylsalicylic acid crystals

Article

Crystal Growth & Design, Vol. 10, No. 7, 2010

3137

Table 4. Elastic Stiffness Coefficients of Acetylsalicylic Acid in GPaa this study atomistic model force field 8

Kim et al.

c11 c22 c33 c44 c55 c66 c12 c13 c23 c15 c25 c35 c46 A λ

9

Ko et al.

experimental

fixed lattice

optimized lattice

DFT

exp

calc

exp

plane wave

RUS

form I

form II

form I

form II

form I

form II

10.76 12.30

11.77

-0.06

10.18 8.14 9.24 3.25 2.13 3.86 5.49 4.50 4.56 -1.73 0.19 -0.64 -0.22

11.55(4) 11.36(3) 10.68(3) 3.67(2) 2.68(2) 4.41(2) 7.73(7) 4.98(4) 5.44(9) -0.49(7) 0.28(9) -0.42(7) -0.08(2)

11.29 12.10 10.67 3.89 2.74 4.52 7.65 4.57 5.67 -0.17 0.60 -0.25 0.32

11.92 11.52 13.63 4.62 3.10 4.47 9.30 7.16 7.74 -2.07 -0.89 -0.81 -1.43

11.93 14.25 14.86 3.95 4.35 2.93 9.75 7.08 6.79 1.57 -0.86 2.20 0.86

13.45 16.30 15.47 4.54 3.95 6.91 12.26 8.59 8.23 -1.23 1.77 -1.62 0.11

13.28 17.90 16.55 3.95 4.55 4.29 12.16 8.19 8.19 1.30 -1.71 2.00 0.37

9.82 7.77 11.71 2.84 3.26 7.41 2.67 4.55 3.97 -1.91 -0.37 -2.57 0.24

12.47 9.70 12.90 5.40 4.72 2.38 7.66 5.20 3.06 1.55 -0.11 3.87 1.76

2.72

7.1 2.13

2.2 2.73

5.4 3.09

2.44

16.2 3.95

2.44

7.8 2.83

1.57

4.57 2.72 3.77 5.88 -0.03 0.64

2.77 4.30

2.68

a The agreement factor A is defined in the text. λ is the smallest eigenvalue of the shear submatrix of the elastic tensor, where the latter is defined as the matrix of the cij with 4 e i, j e 6 in Voigt notation.

Figure 8. Sample holder for resonant ultrasound spectroscopy, with the lower transducer as the ultrasound generator and the upper transducer as the ultrasound detector. are so soft that the sample presumably deforms slightly. Nevertheless, in this setup the crystal in this experiment is close to a freely vibrating body. The resonance spectra were collected at 295.0(5) K in the frequency range from 30 to 800 kHz with a resolution of 0.01 kHz using a network analyzer (Agilent HP 4395A). The sample was mounted in all four possible different orientations, and RUS spectra were collected. In Figure 9, a part of one of the spectra is depicted. The elastic coefficients were obtained using the nonlinear least-squares procedure described above with the aid of a computer program developed by Schreuer18 and are listed in Table 4.

Computational Results The lattice parameters obtained for form I and form II from the COMPASS-based force field geometry optimizations are given in Table 1. The structure of acetylsalicylic acid form I is well described, as the computed lattice parameters agree within ≈3%. The agreement between the COMPASS calculations for form II and the structure suggested by Visheshwar et al.2 is even better. In our calculations form I is less stable by 1.185 kJ/mol than form II.

Figure 9. Detail of a RUS spectrum of acetylsalicylic acid. The complete spectrum covers a frequency range from 30 to 800 kHz.

The computed elastic stiffness coefficients are listed in Table 4, where results are given for fully relaxed structures and for those where the lattice parameters were fixed to the experimental values. There are some discrepancies between the results of the model calculations and the experimental results for form I, but the agreement is satisfactory enough as to imply that the predictions for form II will be semiquantitatively correct. We employed the same measure to estimate if the models are unstable with respect to shear as used by Ouvrard and Price,7 namely the lowest eigenvalue of the “shear submatrix”, i.e. of the 3  3 matrix formed by the coefficients cij with 4 e i,j e 6. The lowest eigenvalues from the experimental data were ≈2.7 GPa. With one exception, the lowest eigenvalue of the shear coefficient submatrix for all models was similar to the experimental values. Only the DFT model of form II gave a lower eigenvalue of 1.6 GPa. Given the agreement between the theoretical and the simulation results, we conclude that form II would not be unstable against shear.

3138

Crystal Growth & Design, Vol. 10, No. 7, 2010

Bauer et al.

Figure 11. Representation surface of the linear thermal expansion of acetylsalicylic acid at 293 K (image created with WinTensor32).

Figure 10. Prediction of the morphology of acetylsalicylic acid crystals from force field based calculations.

The growth morphology can be predicted by calculations of the attachment energies. Attachment energies were calculated for a series of suitable slices drawn from hkl combinations such that the modulus of each index did not exceed 3 and the dhkl distance was not less than 1.3 A˚. The attachment energy is computed as the difference between the lattice energy and the energy of a slice.30 The growth rate is assumed proportional to the attachment energy and a center to face distance is calculated for each face. The morphology is generated using a Wulff plot.31 The predicted morphology for forms I and II is shown in Figure 10. For both forms the faces {100}, {110}, and {011} are among the dominant faces. Additionally, the faces {001} are predicted to exist for form I, while for form II the faces {102} and {111} are predicted to be among the dominant faces. These results for form I are in very good agreement with the observations, as in fact the four predicted faces are the dominant faces of nearly all of the crystals grown from ethanol. It should be possible to tell apart the two polymorphs by their morphologies, if the crystals are growing under formation of the above-mentioned faces, as the angles between these crystal faces are different for the two forms. Discussion The thermal expansion coefficients presented in Table 3 are visualized in Figure 11. There is a significant anisotropy, and the thermal expansion in direction eB2 is about half as large as the effect perpendicular to this direction. The elastic stiffness properties obtained by experiment and atomistic model calculations are given in Table 4, where they are contrasted to those taken from the literature.8,9 The discrepancies between the values obtained by the plain waves/ parallel plate technique and the RUS measurements are generally very small (≈ 0.2 GPa) and only in one case significant (for c22, where the difference is 0.8 GPa). We believe that these minor discrepancies are mainly due to the small deformation of the sample in the RUS measurements. In addition, the

Figure 12. Representation surface of the longitudinal elastic stiffness of acetylsalicylic acid at 295 K. Left: plane wave/parallel plate technique. Right: RUS measurements. Images were created with WinTensor.32

interpretation of RUS spectra of very compressible substances is not straightforward due to overlap of the relatively broad resonance peaks. Specifically, the coefficient c22 can be determined directly from the measurement of a (010) plate with the plane waves/parallel plate technique, and hence, we consider this value to be significantly more reliable than the corresponding RUS-value. The errors associated with the cij obtained from the plane waves/parallel plate technique are due to the leastsquares refinement. For the calculation of the coefficients cij with i ¼ 6 j, more than one data set has to be employed, and that leads to larger absolute errors. However, as all errors are