Elastic Constant Calculations for Molecular Organic Crystals

CRYSTAL GROWTH & DESIGN 2001 VOL. 1, NO. 1 13-27

Articles Elastic Constant Calculations for Molecular Organic Crystals Graeme M. Day, Sarah L. Price,* and Maurice Leslie† Centre for Theoretical and Computational Chemistry, Department of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom Received August 7, 2000

ABSTRACT: Elastic constants of a set of molecular organic crystals have been calculated within the crystal modeling program DMAREL, which was developed to allow the use of highly accurate, anisotropic atom-atom potentials. A set of six molecular crystals (durene, m-dinitrobenzene, the β form of resorcinol, pentaerythritol, urea, and hexamethylenetetramine) has been chosen to sample a range of intermolecular interactions and crystal symmetries. The sensitivity of such calculations to variations in empirical repulsion-dispersion parameters and the electrostatic model is determined and discussed relative to the other errors in the theoretical model and typical experimental uncertainties. We find that model potentials whose functional form has been simplified often reproduce crystal structures and lattice energies adequately but perform poorly when used to calculate elastic constants. This is because more theoretically justified potentials are required to satisfactorily model the curvature at the equilibrium geometries. The rigid-molecule approximation can result in exaggerated elastic constants, and the neglect of thermal effects also leads to significant overestimation of the stiffness constants. Introduction The mechanical properties of crystalline organic products determine many aspects of their commercial production. Thus, the elastic constants of molecular organic crystals are fundamental to the tableting properties of pharmaceuticals1-3 and the assessment of novel (e.g. nonlinear optically active) materials for industrial use.4-6 In addition to the pragmatic interests in predicting elastic properties of industrially important materials, elastic constants are also needed for the interpretation of accurate X-ray-determined electron densities.7-11 Unlike metallic, covalent, and ionic crystals, there is a scarcity of experimentally determined elastic constants for molecular organic crystals.12 This is a consequence of difficulties in growing crystals of adequate size and quality, the predominance of low-symmetry space groups of molecular crystals (which increases the number of independent elastic constants to be determined), and measurement problems resulting from the relative softness of intermolecular interactions. There is, there* To whom correspondence should be addressed. Tel: 020 7679 4622. Fax: 020 7679 7463. E-mail: [email protected]. † Current address: CCLRC Daresbury Laboratory, Warrington WA4 4AD, United Kingdom.

fore, a need to develop and assess computational methods for predicting elastic constants of such crystals. In this paper, we examine the sensitivity of a potentially routine method of calculating elastic constants to the approximations involved and compare the modeling and experimental errors. Elastic properties are largely determined by the forces between molecules, and hence the choice of intermolecular potential is a key factor in the accuracy of calculations. Indeed, while the ability of a given potential to model the solid state is often judged on how well structural parameters and lattice energies are reproduced, the computation of elastic properties of crystals has sometimes been used in the refinement of model intermolecular potentials13,14 and is useful for evaluating existing models. We have calculated the elastic constants using a range of model potentials for a test set of molecules (Table 1), chosen to cover a range of intermolecular interactions. Thus, the sensitivity of elastic constant calculations to variations in the model potential, and other approximations, can be assessed for different types of molecular crystals.

10.1021/cg0055070 CCC: $20.00 © 2001 American Chemical Society Published on Web 11/01/2000

14

Crystal Growth & Design, Vol. 1, No. 1, 2001

Day et al.

Table 1. Molecules Studied, Experimental Crystal Structures, and Elastic Constant Details

a

RT ) room temperature.

Methods and Computational Details The six representative molecules (Table 1) were chosen within the limitations of the most recent compilation of experimental elastic constants.15 Durene (1,2,4,5-tetramethylbenzene) was chosen as an example of a weak, van der Waals crystal, where the electrostatic contribution to the lattice energy is relatively small. m-Dinitrobenzene was included as a molecule with polar functional groups, which increase the importance of the electrostatic contribution. To examine the influence of models for highly directional hydrogen-bonding interactions, we have chosen three molecular crystals exhibiting different hydrogen bonding patternssresorcinol (the β polymorph), pentaerythritol, and urea. The sixth crystal in our set is the very high symmetry body-centered cubic hexamethylenetetramine crystal. The near-spherical shape of this molecule leads to similar packing to ionic or atomic crystals. Although the molecule has polar bonds, steric interactions between methylene groups dominate the crystal packing.

For small organic molecules it is often assumed that the molecular and lattice vibrations are uncoupled. This is justified by the difference in magnitude between inter- and intramolecular forces. When the two types of forces are sufficiently separated, the molecule can be treated as a rigid building block of the crystal lattice. In all of our calculations, we have used the rigid-body approximation and fixed the molecular geometry from the experimentally determined crystal structure. The starting points for the energy minimizations were chosen as the experimentally determined crystal structures and molecular geometries. Because of the difficulty in accurately locating hydrogen atoms from X-ray data, bond lengths to hydrogen atoms in crystal structures determined by X-ray diffraction (durene and m-dinitrobenzene) were adjusted to standard lengths (C-H ) 1.08 Å, N-H ) 1.01 Å, O-H ) 1.02 Å).16 Model Potentials. In addition to the validity of theoretical approximations, the reliability of molecular simulations de-

Elastic Constants for Molecular Organic Crystals

Crystal Growth & Design, Vol. 1, No. 1, 2001 15

Table 2. Buckingham (exp-6) Potential Parameters, Urep-disp(Rik) ) AιK exp(-BιK Rik) - CιKR-6 ik potentiala FIT42

atom pair (W84)50

W9943

GF9351

C‚‚‚C N‚‚‚N (W84)40 O‚‚‚O (W81)41 Hc‚‚‚Hc(W84)50 Hp‚‚‚Hp (FIT)42 C(4)‚‚‚C(4) C(3)‚‚‚C(3) H‚‚‚H C‚‚‚C N‚‚‚N O‚‚‚O H‚‚‚H C‚‚‚N C‚‚‚O C‚‚‚H N‚‚‚O N‚‚‚H (C-)H‚‚‚O (nonpolar) (O-)H‚‚‚O (alcohol) (N-)H‚‚‚O (amine)

A (kJ mol-1)

B (Å-1)

C (Å6 kJ mol-1)

Re (Å)b

(kJ mol-1)b

369 743 254 529 230 064 11 971 5 030 123 404 238 577 13 084 226 145 365 263 195 309 24 158 491 494 393 087 120 792 268 571 228 279 295 432 18 868 777 15 095 064

3.60 3.78 3.96 3.74 4.66 3.60 3.60 3.56 3.47 3.65 3.74 4.01 3.86 3.74 4.10 3.86 4.52 4.82 7.78 7.78

2439.8 1378.4 1123.6 136.4 21.5 1008.0 1642.7 276.0 2418.4 2891.1 1334.7 109.2 2790.7 2681.9 472.8 1523.0 502.1 439.3 1246.8 995.8

3.90 3.66 3.41 3.31 2.48 3.78 3.87 3.26 3.89 3.70 3.61 3.36 3.49 3.61 3.29 3.50 2.98 2.80 1.80 1.80

0.398 0.324 0.400 0.053 0.044 0.194 0.277 0.111 0.387 0.629 0.336 0.042 0.851 0.674 0.205 0.464 0.394 0.505 21.03 16.78

a The W81, W84, and W99 parameters were fitted in combination with an atomic point charge model, the additional FIT parameters with a DMA model, and the GF93 parameters with no electrostatic model. In W99, C(3) and C(4) represent three- and four-coordinated carbon atoms. b Re is the exp-6 potential minimum, and is the potential well depth.

pends on the accuracy with which the intermolecular model potential describes the interaction energies between molecules. Unlike the extreme accuracy with which interaction energies are known between inert-gas atoms36 and the progress toward very accurate, theoretically justified potentials for small polyatomic molecules,37 model potentials describing interactions between typical organic molecules still have a major empirical component38 and are based on the pairwise additive atom-atom potential method, which has been the basis of modeling many crystal properties.39 The atom-atom exchange-repulsion and dispersion contributions to the intermolecular energy are often parametrized simultaneously to define an empirical potential, Urep-disp, often of the exp-6 potential form:

Urep-disp,MN(R,Ω) )

∑

Urep-disp,ik(Rik) )

i∈M,k∈N

∑

Aικ exp(-BικRik) -

i∈M,k∈N

Cικ Rik

6

(1)

where i and k are the constituent atoms of molecules M and N, respectively, and the subscripts ι and κ identify the atom types of atoms i and k. The atom types are chosen according to element and sometimes further subdivided according to bonding environment. We have investigated three of the many available sets of general, transferable exp-6 parameters developed for modeling crystalline organic molecules (Table 2). The W84 potential developed by Williams and Cox40 was parametrized, in conjunction with atomic and off-nuclear point charges, to a large set of crystal structures of oxohydrocarbons41 and azahydrocarbons.40 To decrease the number of parameters needed to describe the four atom types (C, Hc, O, N), a commonly accepted set of combining rules were used to generate the heteroatomic potentials from the homoatomic parameters:

Aικ ) xAιιAκκ, Bικ )

BιιBκκ , Cικ ) xCιιCκκ 2

(2)

The W84 set of parameters was later supplemented by exp-6 parameters for polar hydrogens (Hp) by fitting to a large set of hydrogen-bonded organic molecules using a distributed multipole electrostatic model.42 We refer to this expanded potential set as FIT. Williams later parametrized an improved set of hydrocarbon parameters, W99,43 which distinguish between aliphatic and aromatic carbons by fitting to a larger

collection of hydrocarbon crystal structures and lattice energies, in conjunction with a multicenter point-charge electrostatic model. In this work, the lattice energy contributions from the exp-6 terms were summed to a 15 Å cutoff. These empirical potentials were parametrized to represent all intermolecular interactions except for an explicit electrostatic model. The electrostatic models are designed to be sufficiently accurate that we assume that the exp-6 parameters are independent of the electrostatic model. We consider various electrostatic models, all derived from an MP2(full)/6-31G(d,p) wave function of the isolated molecule, calculated within the Gaussian98 suite of electronic structure programs.44 Secondorder Møller-Plesset perturbation theory45 describes electron correlation well enough to avoid the exaggerated charge separation that is characteristic of self-consistent-field wave functions. The polarized split valence basis set provides sufficient angular and radial flexibility in the basis functions to adequately describe the charge density in polar organic molecules. The simplest, least computationally intensive electrostatic model involves assigning atomic point charges. There are many methods of obtaining atomic charges, ranging from Mulliken46 charges to electrostatic potential (ESP) derived charges, which are optimized to reproduce the molecular electrostatic potential on grids outside the molecule. The usefulness of different schemes for assigning atomic charges has been examined elsewhere.47,48 We have chosen the CHelpG algorithm of Breneman and Wiberg,49 with charges constrained to reproduce the molecular dipole moment and preserve any molecular symmetry. These charges were calculated within the Gaussian98 suite of electronic structure programs.44 In some cases this model was contrasted with Mulliken charges to test the sensitivity to the method used to obtain atomic charges. While the point-charge model is a good starting approximation, some features that arise from the rearrangement of valence electron density upon bonding, such as localized lone pairs and π-electron density, cannot be adequately described using an isotropic description of atoms. A rigorous approach to describe the electrostatic interactions involves the expansion of the electron density in terms of atom-centered multipoles: e.g., placing a charge, dipole, quadrupole, etc. at each atomic site. Several methods have been proposed for generating such distributed multipoles from the molecular wave function.52-56 We used the distributed multipole analysis (DMA),55,56 which optimizes the convergence of the multipole series, as implemented in the electronic structure package CADPAC,57 to generate atom-centered distributed multipoles, up to hexadecapole, for all of the molecules in our test set. This anisotropic

16


Day et al.

electrostatic model results in a more complicated representation of the atom-atom interaction energy, which now depends on the relative orientation of the atoms as well as the interatomic separation, and produces noncentral forces and torquesbetweentheatoms.Themechanicsofsuchinteractions58-60 are programmed into ORIENT60 for modeling molecular clusters and DMAREL59 for molecular crystals. The increased accuracy in the electrostatic model has been shown to be vital in modeling the geometries of van der Waals clusters61-63 and the structures of many molecular crystals.64 This is especially true where hydrogen bonding is involved, as the electrostatic interaction of the polar hydrogen atom with localized lone pairs usually determines the optimal orientation of such interactions. The electrostatic contributions to the lattice energy were evaluated by the Ewald summation technique65 for chargecharge, charge-dipole, and dipole-dipole terms. All higher order electrostatic interactions, up to R-5 ik in the multipole expansion, were summed to a 15 Å cutoff between the centers of mass of the molecules. Thus, the DMA electrostatic model gives quite an accurate estimate of the electrostatic contribution to the lattice energy. However, the polarization and electrostatic penetration terms are not explicitly represented in the model. The problems and expense involved in obtaining and using explicit electrostatic models led Filippini and Gavezzotti66,67 (GF93) to develop an exp-6-only potential model by fitting to the crystal structures and sublimation energies of a large data set of simple monofunctional organic molecules. To allow the exp-6 parameters to absorb the electrostatic effects as far as possible, Filippini and Gavezzotti abandoned the use of combining rules for heteroatomic interactions. As is evident in Table 2, the resulting parametrization leads to interatomic potentials with vastly different well depths, , and equilibrium distances, Re, from the empirical potentials derived with an explicit electrostatic model. Elastic Constant Calculations. Several methodologies have been used to calculate elastic constants from the second derivatives of the potential energy surface at the classical, 0 K structure. These methods include full lattice dynamic calculations, which involve extracting wave velocities from calculated phonon dispersion curves and relating these to the elastic constants,68-72 the numerical application of axial and shear strains,73 and analytic expressions directly relating second derivatives of the potential energy hypersurface to the elastic constant matrix.74-78 We have incorporated the analytic method for rigid molecular crystals modeled by atom-atom interactions into the crystal structure modeling program DMAREL,59 thus enabling the use of highly accurate, anisotropic electrostatic models in elastic constant calculations. Using this methodology,74,75 the increase in energy density, ∆Ψ, due to a homogeneous deformation (a deformation in which the lattice structure of the crystal is retained) is expanded in a Taylor series and truncated at harmonic terms. Forces, strains, and torques vanish at the equilibrium crystal structure, and the second derivatives of the lattice energy are related to the internal and external strain contributions to ∆Ψ. The increase in energy density is then related to the stiffness constants by

∆Ψ ) 1/2

∑C

ij i j

(3)

i,j

where i are elements of the strain matrix of elastic theory.79 As the lattice energy is calculated directly from the potential energy surface, zero-point and thermal effects are ignored. Thus, the resulting elastic constants correspond to the classical, 0 K state. Starting from the experimental crystal structure and using a modified Newton-Raphson procedure, analytical forces, rigid molecule torques, and strain derivatives of the unit cell59 were used to locate a stationary point on the potential energy surface defined by the chosen model potential. The secondderivative matrix at the stationary point was calculated

numerically from the analytical gradients. The minimization path was initially constrained to the observed symmetry of the crystal, and the eigenvalues of the Hessian were then used to characterize the stress-free structure as a minimum or saddle point on the potential energy surface. If a true minimum could not be found in the experimental space group, then the symmetry was relaxed to a subgroup, chosen according to the negative eigenvalues, until an energy minimum was found. Thus, the calculated elastic constants correspond to the stress-free crystal structure at a minimum point on the potential energy surface defined by the atom-atom potential.

Results Tables 3-8 present the equilibrium lattice parameters and calculated elastic constants of the six molecular crystals in Table 1, calculated for a range of model potentials. Experimental structures and elastic constants are presented for comparison. In discussing the calculated elastic constants, we focus on the errors in the diagonal elements of the elastic constant matrix (Cjj, j ) 1-6), describing the stiffness to uniaxial compression and shear. Because of the manner in which elastic constants are determined from raw experimental data, the experimental uncertainty is lowest in these values.69,80 We report our calculated off-diagonal elements (Cij, i * j), which describe biaxial compression and distortion of the crystal, as they are also important for describing mechanical properties of the crystal. However, large experimental errors often preclude meaningful comparisons. Durene, C6H2(CH3)4. Because of the validity of the rigid-body approximation and the simplicity of the model potential, there have been many calculations of the elastic constants of aromatic hydrocarbons, such as benzene,73,77,81 naphthalene,77,82 and anthracene.82 Early calculations ignored the electrostatic contribution to intermolecular interactions, which, for unsubstituted hydrocarbons, is a minor contribution relative to the van der Waals (exchange-repulsion and dispersion) energies. In fact, using the W99 + DMA potential model, we estimate that electrostatic interactions account for less than 13% of the total lattice energy of durene. Durene crystallizes in the monoclinic space group P21/c (Z ) 2) with half a molecule in the asymmetric unit. The molecules pack in a herringbone motif with close contacts between methyl hydrogens and the face of the aromatic ring of a neighboring molecule (Figure 1). These contacts, between nearly perpendicular molecules in the a and b directions, are stabilized by dispersive CH3‚‚‚π interactions. There are layers stacked in the c direction, where the only close contacts are between methyl and aromatic hydrogen atoms. All of the model potentials reproduce the structure of the durene crystal moderately well (Table 3a), with an rms percent error in the lattice parameters of less than 5%. There is very little difference in the final structures using the GF93 and exp-6 + ESP potentials. While enhancement of the electrostatic model (exp-6 + DMA) improves the description of methyl-aromatic ring interactions slightly, as reflected in the a and b lattice parameters, the methyl-methyl interactions in the c direction are modeled less well, allowing the structure to collapse by more than 7% in this direction. The measured stiffness constants show that the crystal is nearly isotropic to compression (i.e. C11 ≈ C22



Figure 1. Crystal structure of durene, showing crystallographic and optical axes. Figures have been prepared using the Cerius2 software.83 Table 3. Experimental and Calculated Structures, Lattice Energies, and Elastic Constants of Durene (a) Structures and Lattice Energiesa a (Å) exptl17 W6718 GF93 W84 + ESP W99 + ESP W84 + DMA W99 + DMA

b (Å)

11.59

5.74

11.78 (+1.61) 12.10 (+4.37) 12.09 (+4.29) 11.48 (-0.98) 11.48 (-0.69)

5.43 (-5.49) 5.55 (-3.25) 5.63 (-1.92) 5.90 (+2.75) 5.96 (+3.83)

c (Å)

V (Å3)

β (deg)

rms error (%)b

-Φlatt (kJ mol-1)c 71.7 ( 0.3

7.04 112.8 431.8 no structural or energy details given 6.97 (-1.03) 114.51 (+1.52) 405.0 (-6.19) 6.93 (-1.60) 114.11 (+1.16) 424.7 (-1.62) 6.96 (-1.07) 113.77 (+0.86) 433.7 (+0.46) 6.51 (-7.56) 105.93 (-6.09) 423.6 (-1.89) 6.53 (-7.31) 105.21 (-6.73) 431.9 (+0.04)

3.36 3.28 2.79 4.68 4.78

73.72 71.96 63.38 78.24 69.47

(b) Elastic Constants (GPa) exptl18 W6718 GF93 W84 + ESP W99 + ESP W84 + DMA W99 + DMA

exptl18 W6718 GF93 W84 + ESP W99 + ESP W84 + DMA W99 + DMA

C11

C22

C33

C44

C55

C66

9.08 9.2 12.43 10.82 9.15 12.45 10.68

10.00 7.3 14.27 11.66 9.75 12.59 10.71

10.03 6.6 13.72 12.26 10.54 16.77 14.30

1.84 2.0 2.57 2.17 2.03 4.23 3.75

2.17 4.3 5.09 4.92 4.06 3.20 2.73

7.27 10.7 13.38 11.44 9.68 12.19 10.44

C12

C13

C23

C15

C25

C35

C46

7.75 4.5 10.57 8.95 7.33 9.44 7.67

2.57 2.4 8.60 8.20 7.19 7.08 6.08

3.19 1.0 6.23 6.09 5.28 8.05 7.31

0.11 -0.8 0.45 0.44 0.56 -0.10 0.19

1.04 0.4 -2.41 -2.21 -1.85 -1.70 -1.34

-0.13 -0.0 2.41 2.53 2.15 2.75 2.18

0.13 0.4 -3.08 -2.73 -2.31 -3.11 -2.55

a Percentage errors given in parentheses. b RMS of percentage errors in a, b, and c. c The experimental lattice energy, Φ latt, is approximated by -∆H°sub.84

≈ C33), while the diagonal shear constants (C44, C55) indicate a weakness to shearing perpendicular to the c axis. All of the parameter sets reproduce these qualitative aspects of the elastic stiffness constants fairly well.

Sanquer and Ecolivet’s calculations used the exp-6 potential parametrized to hydrocarbons by Williams (W67)85 with no explicit Coulombic term. Their results agree with the experimentally determined elastic con-

18


Day et al.

Figure 2. Crystal structure of m-dinitrobenzene, showing crystallographic and optical axes.

stants reasonably well, but the stiffness to compression in the b and c directions is significantly underestimated. All other potentials overestimate the stiffness to compression. The rms percent errors in the diagonal elements of the stiffness matrix are 72.2% using the GF93 potential, 48.0% (Sanquer and Ecolivet,18 W67), 58.9% (W84), and 38.3% (W99) with ESP point-charge models and 71.0% (W84) and 50.8% (W99) using DMA electrostatics. Sanquer and Ecolivet18 suggested that thermal motion of the methyl hydrogens may account for much of the error in the calculated elastic constants. Among our calculations, the change in potential that results in the most noticeable improvement is the distinction between three- and four-coordinated carbon atoms in the W99 exp-6 parametrization, independent of the electrostatic model. The quality of the results using a DMA electrostatic model is worse than with ESP charges, but this can be attributed to the problems in reproducing the crystal structure with the exp-6 + DMA potentials. m-Dinitrobenzene, C6H4(NO2)2. m-Dinitrobenzene, which crystallizes in the orthorhombic space group Pbn21 (Z ) 4), is a molecule of size and shape similar to those of durene. Like durene, m-dinitrobenzene forms a herringbone motif and the dominant contacts in the bc plane are the edge-to-face contacts with the nitro groups pointed into the π-charge density of the neighboring molecule (Figure 2). We estimate that the repulsion-dispersion contributions to the lattice energies of the two crystals are almost identical. However, due to the local polarity of the nitro substituents, the electrostatic contribution to the energy of the m-dinitrobenzene crystal is approximately double the contribution to the lattice energy of durene. The use of an explicit electrostatic model, even in its simplest point-charge form, reduces the rms percent

error in the lattice constants by half, while also bringing the lattice energy into much better agreement with the experimental heat of sublimation (Table 4a). The distributed multipole model, with the same exp-6 parameters, gives a further improvement of the calculated lattice parameters, resulting in an rms percent error of less than 1%. Because of the similarity in crystal structure, the stiffness in the x and y directions (C11, C22, and C12) is similar to that in durene (Table 4b). However, the interaction of the local dipoles of the nitro groups, as well as close N-O‚‚‚H-C contacts, results in a uniaxial stiffness constant, C33, double that of durene in the z direction, orthogonal to the edge-to-face contacts. Overall, the GF93 potential greatly overestimates the stiffness of the crystal. Both potentials which include Coulombic terms give improved agreement with the experimentally determined stiffness constants, with a particular improvement in the shear constant, C55, which shows a greater stability toward xz shearing, in agreement with experiment. The error in the diagonal stiffness constants is greatly reduced by treating electrostatic interactions explicitly (rms errors of 38.2% using W84 + ESP, 37.3% using W84 + DMA, compared with 94.8% from GF93 calculations). Replacing atomic charges by DMA atom-centered multipoles results in only modest improvement. Resorcinol (β Form), C6H4(OH)2. The β polymorph of resorcinol (m-dihydroxybenzene) forms a crystal structure similar to that of durene and m-dinitrobenzene, in the orthorhombic space group Pna21 (Figure 3). However, resorcinol molecules have strong hydrogen bond donor and acceptor sites, leading to a network of intramolecular hydrogen bonding. While there are weak hydrogen-bonding interactions in the a direction, between the buckled molecular sheets, the major components of these interactions are oriented in the b and



Table 4. Experimental and Calculated Structures, Lattice Energies, and Elastic Constants of m-Dinitrobenzene (a) Structures and Lattice Energiesa exptl19 GF93 W84 + ESP W84 + DMA

a (Å)

b (Å)

c (Å)

V (Å3)

rms error (%)b

-Φlatt (kJ mol-1)c

13.26 12.91 (-2.63) 13.12 (-1.05) 13.14 (-0.87)

14.05 13.97 (-0.54) 13.80 (-1.73) 13.91 (-0.98)

3.81 3.65 (-4.08) 3.85 (+1.06) 3.80 (-0.13)

708.8 658.5 (-7.10) 696.5 (-1.73) 694.8 (-1.97)

2.82 1.32 0.76

81.2 ( 1.7 105.38 82.97 84.38

(b) Elastic Constants (GPa) exptl21 exptl20 GF93 W84 + ESP W84 + DMA

C11

C22

C33

C44

C55

C66

C12

C13

C23

10.70 14.8 19.88 13.24 13.90

11.30 12.7 26.46 16.87 16.71

20.27 18.1 37.80 28.22 28.44

4.37

2.04

5.31

6.30

1.95

3.19

9.10 5.98 6.32

0.35 1.10 1.51

7.96 6.75 6.88

15.22 10.75 10.45

10.55 6.03 5.83

14.41 7.75 7.11


Figure 3. Crystal structure of β-resorcinol, showing crystallographic and optical axes. Hydrogen bonds are shown with fine green dashed lines.

c directions. From our calculations, electrostatic interactions account for approximately two-thirds of the lattice energy. The relaxed structures (Table 5a) demonstrate that the FIT + Mulliken model gives the worst overall reproduction of the experimental crystal structure, particularly in the a and c directions. Of the models including explicit electrostatics, the errors in the equilibrium lattice parameters decrease dramatically as the form of the electrostatic model is improved. Althoughthe distorted, equilibrium structures are quite different, the FIT + ESP potential has approximately the same 7% rms error in the lattice constants as the GF93 relaxed structure, but the better DMA electrostatic model leads to an equilibrium geometry that varies only slightly from the experimental structure (rms error 1.21%). A similar trend is observed in the calculated elastic constants (Table 5b). As with durene and m-dinitrobenzene, the GF93 potential results in unacceptable

errors, with the stiffness constants being generally overestimated. The elastic constants calculated using Mulliken charges are also in poor agreement with experiment, but the errors are not so systematicssome constants are predicted to be too soft, while others are too stiff. The errors are decreased significantly with the use of ESP electrostatics and further improved when the atomic charges are replaced by atom-centered multipole expansions. Pentaerythritol, C(CH2OH)4. Pentaerythritol crystallizes in a body-centered tetragonal space group, (I4 h, Z ) 2) with hydrogen bonding extending in the ab plane (Figure 4). Each molecule accepts and donates a hydrogen bond to each of its four nearest neighbors in the molecular sheets. Hence, each molecule is held in place by a total of eight hydrogen bonds. The strengths of interactions in the ab plane are, thus, much stronger than in the c direction, perpendicular to the molecular sheets, in which the closest contacts between the molecules are CH2‚‚‚H2C steric interactions.

20


Day et al.

Table 5. Experimental and Calculated Structures, Lattice Energies, and Elastic Constants of β-Resorcinol (a) Structures and Lattice Energiesa exptl22 GF93 FIT + Mulliken FIT + ESP FIT + DMA

a (Å)

b (Å)

c (Å)

V (Å3)

rms error (%)b

-Φlatt (kJ mol-1)c

7.81 6.85 (-12.27) 9.37 (+19.97) 8.62 (+10.37) 7.78 (-0.37)

12.62 12.71 (+0.74) 12.94 (+2.59) 12.93 (+2.46) 12.87 (+1.98)

5.43 5.74 (+5.75) 4.54 (-16.35) 4.99 (-7.99) 5.46 (+0.57)

534.7 499.8 (-6.54) 550.5 (+2.95) 556.4 (+4.04) 546.5 (+2.19)

7.84 14.98 7.69 1.21

87.5 ( 0.5 108.94 80.09 86.26 93.76

(b) Elastic Constants (GPa) exptl23 GF93 FIT + Mulliken FIT + ESP FIT + DMA

C11

C22

C33

C44

C55

C66

C12

C13

C23

8.6 17.91 15.73 12.03 9.93

28.8 42.11 33.11 34.66 35.22

19.5 25.86 9.66 11.47 14.98

3.26 21.27 4.77 5.96 9.83

4.35 8.21 5.99 7.52 7.79

4.00 2.90 4.11 2.73 3.54

9.5 13.32 11.35 12.17 9.63

7.5 18.03 10.84 11.16 10.60

19.1 19.89 9.04 8.84 11.00

a Percentage errors given in parentheses. b RMS of percentage errors in a, b, and c. c The experimental lattice energy, Φ latt, is approximated by -∆H° sub for the R polymorph.87

Figure 4. Crystal structure of pentaerythritol, showing crystallographic and optical axes. Hydrogen bonds are shown with fine green dashed lines.

The a and b lattice parameters are reproduced very well by all three sets of potentials, allowing the unit cell to expand by less than 1.5% in these equivalent directions (Table 6a). The dominant interactions in this crystal plane are the hydrogen bonds, whose equilibrium (O‚‚‚HP) lengths are only 2.9% too long for FIT + DMA, compared to 4.7% for FIT + ESP and 3.1% for GF93. The length of the unit cell in the c direction is mainly dictated by interlayer CH2‚‚‚H2C van der Waals interactions. All of the potentials allow this lattice parameter to contract when the stresses and forces are relaxed, and this effect is most pronounced in the final GF93 structure, where the distance between molecular sheets has decreased by more than 10% relative to the experimental structure. The FIT + ESP and FIT + DMA potentials reproduce the cell length in the c direction satisfactorily, with errors less than 4%. The two more recent experimental determinations of the stiffness matrix are in good accord, while the values

presented by Srivastava and Chakraborty,25,26 as determined from the intensity of diffuse X-ray reflections, are in gross disagreement. The three general, transferable model potentials severely overestimate the stiffness to uniaxial stress in the ab plane, described by C11, with the GF93 potential predicting a value almost an order of magnitude higher than the experimental results and FIT + DMA giving a 3-fold exaggeration (Table 6b). The magnitude of this error is much greater than for any of the other molecular crystals considered in this study. In contrast, the calculated values using a specifically developed potential89 are in excellent agreement with experiment. Morse potentials and empirically fitted atomic charges were used to describe the inter- and intramolecular interactions in the pentaerythritol crystal, thus allowing for molecular flexibility in their calculations. This 34-parameter model was successfully parametrized to reproduce the infrared and Raman vibrational spectrum of the pentaerythritol crystal to



Table 6. Experimental and Calculated Structures, Lattice Energies, and Elastic Constants of Pentaerythritol (a) Structures and Lattice Energiesa exptl24 GF93 FIT + ESP FIT + DMA

a ) b (Å)

c (Å)

V (Å3)

rms error (%)b

-Φlatt (kJ mol-1)c

6.08 6.14 (+0.97) 6.17 (+1.47) 6.14 (+0.94)

8.74 7.83 (-10.41) 8.42 (-3.71) 8.44 (-3.50)

323.2 295.2 (-8.66) 320.4 (-0.85) 317.7 (-1.69)

6.06 2.32 2.16

163 162.88 133.41 142.42

(b) Elastic Constants (GPa)a exptl25,26 exptl27 exptl28 calcd89 GF93 FIT + ESP FIT + DMA

C11

C33

C44

C66

C12

C13

C16

2C66 + C12

6.1 ((4) 40.5 ((1) 39.4 ((0.25) 39.95 317.14 110.82 117.88

8.0 ((4) 13.9 ((1) 13.3 ((1.5) 13.67 31.06 16.98 17.43

3.5 ((4) 2.74 ((1) 4.4 ((7) 2.58 15.57 8.00 9.56

4.6 ((4) 2.52 ((1)

-2.50 ((6) 26.6 ((8)

0.50 ((6) 10.5 ((2.5) 7.6 ((4) 9.13 9.23 6.88 7.91

-0.39 ((4) 3.13 ((76) 0.0 ( 0.5 1.97 -2.71 1.58 1.22

6.7 ((8) 31.6 ((7) 5.13 ((0.6) 29.8 32.07 16.93 27.41

2.49 9.38 8.92 13.92

24.82 0.93 -0.91 -0.43


within 4-8% in intramolecular (CH2 scissoring and OH torsional) modes and 2% in all other frequencies.90 Considering the close relationship between elastic constants and zone center vibrations, the success of the stiffness constant calculations is not surprising when parameters have been optimized to reproduce experimental frequencies. These results suggest that intramolecular flexibility plays an important role in the elastic behavior of pentaerythritol. In fact, the Raman,91 infrared,92 and (k ) 0) calculated90 vibrational spectra of the pentaerythritol crystal all have several low-energy (165-425 cm-1) bands that have been assigned to C-C-O motions. These low-energy intramolecular modes cannot be treated as uncoupled from the highest energy lattice vibrations at approximately 75 cm-1. Hence, we cannot assume that the rigid-body approximation is appropriate for this crystal. The large errors in our elastic constant calculations are consistent with the inapplicability of the rigid-body approximation for modeling pentaerythritol. Allowing angle bending and C-C-O-H torsional distortion in the simulations would increase the compliance of the crystal in the a and b directions, as well as softening the crystal to shearing forces, which are also badly modeled in the rigid-molecule approximation (C44 and C66). The only stiffness constant that is relatively unaffected by C-C-O flexibility is the C33 uniaxial constant, describing stiffness to changes in the spacing between molecular sheets. While the value of this stiffness constant is overestimated by the GF93 potential, the FIT + ESP and FIT + DMA model potentials both give only a modest overestimate. Urea, CO(NH2)2. While the crystal structures of resorcinol and pentaerythritol are characterized by O-H‚‚‚O hydrogen bonds that extend mainly in two dimensions, urea crystallizes such that it forms a strong three-dimensional array of N-H‚‚‚O hydrogen bonds. The polar molecules align so that there are, in fact, four hydrogen bonds to each carbonyl oxygen, resulting in a particularly stiff molecular organic crystal. The urea molecule (Table 1) adopts C2v symmetry in the tetragonal (P4 h 21m) crystal structure (Figure 5). There have been several previous computational studies of the elastic and vibrational properties of urea.73,93,94 Pavlides et al.73 and George et al.93 used the same intermolecular potential in their calculationssa Len-

nard-Jones 12-6 plus point-charge model parametrized to several amide crystal structures.95 George et al. supplemented this intermolecular potential by empirical functions describing harmonic bond stretching, angle bending, and torsions. The relaxed structures (Table 7a) and calculated stiffness constants (Table 7b) are in excellent agreement, which suggests that the rigidmolecule approximation does not introduce serious errors. Lattice energy minimization has problems in reproducing the tetragonal structure of urea. The equilibrium structure reported by Pavlides and co-workers shows slight deviations from the observed tetragonal symmetry73 and, while stationary points in tetragonal symmetry were located on all four of our potential energy surfaces (GF93, FIT + Mulliken, FIT + ESP, and FIT + DMA), analysis of the Hessian matrix characterized these as transition structures between two symmetry-equivalent minima. The imaginary frequency corresponds to rotation of the urea molecule around its principal axis, breaking the equivalence between the a and b directions and resulting in distortion to an orthorhombic subgroup of the experimental symmetry. Although the geometrical effect is significant, the resulting lowering of lattice energy is negligible (e.g. 1.4 kJ/mol on the FIT + DMA potential surface), indicating a very flat potential energy. As further evidence for this observation, the stiffness matrix given by George et al.93 has one eigenvalue numerically equal to zero (C11 ) C12), corresponding to a zero gradient along the path describing the tetragonal to orthorhombic distortion. These observations are consistent with the soft acoustic mode in the [110] direction near the Brillouin zone center of urea, calculated by Parlinski and Chapuis94 using the Dreiding inter- and intramolecular force field.96 This consistency over so many potential models suggests that the problem is not in the parametrization of the atomatom model potential. Hartree-Fock calculations have estimated that the induction energy accounts for approximately 9% of the lattice energy,97 suggesting limitations of the pairwise additive approximation for crystalline urea. (Although the structural and vibrational properties of crystalline urea have been studied using periodic Hartree-Fock and density functional theory,97-100 these calculations have assumed the experimental symmetry.) Alterna-

22


Day et al.

Figure 5. Crystal structure of urea, showing crystallographic and optical axes. Hydrogen bonds are shown with fine green dashed lines. Table 7. Experimental and Calculated Structures, Lattice Energies, and Elastic Constants of Urea (a) Structures and Lattice Energiesa exptl29 calcd93 calcd73 GF93 FIT + Mulliken FIT + ESP FIT + DMA

a (Å)

b (Å)

c (Å)

V (Å3)

rms error (%)b

5.56 5.55 (-0.18) 5.546 (-0.18) 5.31 (-4.62) 5.69 (+2.17) 5.18 (-6.88) 5.59 (+0.43)

5.56 5.55 (-0.18) 5.551 (-0.18) 4.91 (-11.85) 4.99 (-10.36) 5.39 (-3.12) 4.97 (-10.73)

4.68 4.70 (+0.4) 4.79 (+2.35) 4.48 (-4.41) 4.86 (+3.67) 4.78 (+2.15) 4.73 (+1.00)

145.1 144.8 (-0.20) 147.5 (+1.65) 116.6 (-19.63) 137.7 (-5.05) 133.7 (-7.84) 131.3 (-9.46)

0.29 1.36 7.77 6.47 4.53 6.23

-Φlatt (kJ mol-1)c 98.6 112.56 63.25 86.11 85.80

(b) Elastic Constants (GPa)a exptl30 exptl31 exptl32 calcd93 calcd73 GF93 FIT + Mulliken FIT + ESP FIT + DMA

C11

C22

C33

C44

C55

C66

C12

C13

C23

21.7 ((2) 23.5 ((2) 11.2 17.1 18.19 18.42 39.40 10.37 8.22 14.02

21.7 ((2) 23.5 ((2) 11.2 17.1 18.19 18.42 39.96 10.57 11.92 15.59

53.2 ((1) 51.0 ((2)

6.26 ((1) 6.2 ((2) 5.9 5.1 5.15 5.12 10.02 2.47 7.13 4.49

6.26 ((1) 6.2 ((2) 5.9 5.1 5.15 5.12 13.70 4.48 5.69 7.44

0.45 ((4) 0.50 ((10) 8.44 22.2 23.61 23.65 25.26 8.87 12.00 8.74

8.9 ((21) -0.5 ((10) 10.2 17.1 17.71 17.52 32.33 7.84 7.88 8.43

24.0 ((17) 7.5 ((11)

24.0 ((17) 7.5 ((11)

11.6 12.87 12.58 16.76 5.80 5.25 10.45

11.6 12.87 12.58 12.76 3.85 10.00 4.25

68.7 63.90 61.65 182.81 38.69 50.29 55.47


tively, the calculated symmetry breaking could be a result of the classical, 0 K nature of the simulations. Because of the small energy change corresponding to the orthorhombic distortion, the zero-point motions could average over the symmetrically equivalent, distorted minima, resulting in the experimentally observed tetragonal symmetry. Analyzing X-ray thermal parameters, Buergi101 found considerable deformations of the amine groups, especially at the zone boundaries. Thus, correctly representing the experimental tetragonal structure of urea appears to be a challenge for current modeling methods. Despite the breaking of the symmetry, the main

features of the calculated and experimental stiffness matrices can still be reasonably compared (Table 7b). Because of their mutual agreement, it seems best to compare to the experimental stiffness constants reported by Fischer and Zarembowitch30 and Yoshihara and Bernstein.31 However, Kretchetov et al.32 have determined an incomplete stiffness matrix that is in poor agreement with the other two sets, suggesting that the experimental uncertainty is large. In particular, we note the discrepancies in the elastic constants which relate to stresses in the ab plane, C66 and C12. The two complete sets of experimental elastic constants describe surprisingly little stiffness to xy shearing, but none of



Figure 6. Crystal structure of hexamethylenetetramine, showing crystallographic and optical axes.

the calculated values support this apparent compliance. Instead, the calculations corroborate the value determined by Kretchetov et al.32 The largest axial stiffness constant is observed along the c direction of the crystal. All of the hydrogen bonds have a large component in this direction, parallel to the dipole moment of the urea molecule. The corresponding stiffness constant, C33, is reproduced quite well by most of the model potentials, except the GF93 potential, which predicts a stiffness more than 3 times the observed value, and the FIT + Mulliken potential, which gives a value that is far too soft. Hexamethylenetetramine, C6N4H12. Hexamethylenetetramine is an unusual molecular crystal because of its very high molecular (Td) and crystal symmetry (I4 h 3m, Z ) 2) (Figure 6). According to a 1995 survey,103 only 0.2% of homomolecular organic crystals are found in cubic space groups. Hence, there are only three independent elastic constants, and this high symmetry is one reason hexamethylenetetramine was one of the first molecular organic crystals whose elastic properties were measured.34,35 As with other crystals with no directional intermolecular interactions, the GF93 potential reproduces the structural parameters of the hexamethylenetetramine crystal very well (Table 8a). Furthermore, the minimized structure is only slightly improved when an explicit electrostatic model is introduced, though the calculated lattice energy is brought in better agreement with the experimental heat of sublimation. The GF93 potential overestimates the stiffness of the hexamethylenetetramine crystal significantly, predicting all three stiffness constants to be double the observed values (Table 8b). The potentials with ESP and DMA electrostatic descriptions reproduce the stiffness constants more realistically, approximately halving the errors. The rms percent errors over all three independent elastic constants are 113.4% with the GF93 po-

Table 8. Experimental and Calculated Structures, Lattice Energies, and Elastic Constants of Hexamethylenetetramine (a) Structures and Lattice Energiesa exptl33 GF93 W84 + ESP W84 + DMA

a ) b ) c (Å)

V (Å3)

-Φlatt (kJ mol-1)b

6.927 6.83 (-1.45) 7.01 (+1.25) 7.01 (+1.26)

332.4 318.1 (-4.32) 345.0 (+3.76) 345.1 (+3.82)

74.9 ( 2.9 96.57 78.87 84.08

(b) Elastic Constants (GPa)a exptl35 exptl34 GF93 W84 + ESP W84 + DMA

C11

C12

C44

15 16.43 ((0.3) 41.01 25.71 26.51

3 4.33 ((2.5) 7.68 6.12 6.60

7 5.15 ((1) 10.36 7.81 7.80

a Percentage errors given in parentheses. b The experimental lattice energy, Φlatt, is approximated by -∆H°sub.104

tential, 50.2% with W84 + ESP, and 55.3% using the W84 + DMA model. Discussion We have calculated the elastic constants for six molecular crystals, with a variety of model potentials. Apart from those employing Mulliken charges, the models all give reasonable reproductions of the crystal structures and lattice energies, given the expected limitations of lattice energy minimization. (The neglect of temperature effects may lead to errors comparable with thermal expansion, which is typically a few percent in the cell lengths for organic crystals. The many approximations involved in comparing lattice energies with enthalpies of sublimation mean that differences of around 8-10 kJ/mol are no cause for concern.51) For elastic constants, it is necessary to discuss how the observed variations in the calculated elastic constants with model potential, and their large deviations from the experimental values, are determined by the various approximations in the calculations.

24


Variations in the Model Potential. This study shows that variations in the functional form of the model potential, as well as the parametrization used, have significant effects on the calculated elastic properties of molecular organic crystals. In calculating elastic constants, it appears necessary to treat repulsiondispersion and electrostatic contributions by separate terms in the model potential. The GF93 potential, which absorbs average electrostatic effects into the exp-6 model, was developed specifically to provide a costeffective method of modeling crystal structures and lattice energies, which it does remarkably well. However, the potential fails to reproduce the elastic properties of molecular crystals satisfactorily. Over this set of diverse molecular crystals, the stiffness constants are consistently badly overestimated, particularly in the crystallographic directions where the electrostatic component of the interaction is dominant. These observations indicate that, while the GF93 model reproduces satisfactory well depths and equilibrium distances, the exp-6 functional form cannot correctly describe the curvature of the potential energy surface in the vicinity of the equilibrium structures. This is consistent with the model potential decaying as R-6, whereas electrostatic interactions decay more slowly. Thus, elastic constant calculations clearly demonstrate the problems of using empirically fitted potentials for modeling properties that depend on aspects of the potential energy surface that are not sampled in the parametrization. The improvement in the calculated elastic constants through using a more realistic functional form with an explicit electrostatic model definitely justifies the increase in computational expense and complexity. Molecular crystal structure modeling for all but the saturated hydrocarbons requires a reasonably realistic electrostatic model, as shown by the results using Mulliken atomic charges, where the calculated structures and elastic constants are in poor agreement with experiment. The more rigorous electrostatic models used in the present work are designed to accurately reproduce the electrostatic potential at typical intermolecular contacts. The added complexity of replacing the ESP point charges with the DMA anisotropic atom-atom potential normally improves the agreement of the relaxed structural parameters with the experimental structure. The increased accuracy in the calculated elastic constants is less pronounced but is noticeable when strong, directional interactions are involved, as in the resorcinol, pentaerythritol, and urea crystals. Overall, when a contribution to the model intermolecular potential becomes more theoretically justified, there is an improvement in the elastic constants in the crystal directions where this contribution is significant. This is seen both in the increasing realism of the electrostatic models, from no explicit electrostatic model, through Mulliken and ESP charges to distributed multipoles, for the hydrogen-bonded and polar molecules, and in the improved parametrization of the repulsiondispersion in differentiating between aromatic and aliphatic carbons for the van der Waals crystal, durene. However, even the most realistic model potentials are systematically overestimating the elastic constants of all the molecular crystals considered, to a degree that

Day et al.

cannot be due only to the remaining approximations in the potential. Therefore, some assessment of the likely errors resulting from other approximations in the theory and comparison with experiment is necessary for interpreting these results. Rigid-Body Approximation. The rigid-body approximation led to a gross overestimate of the elastic constants of pentaerythritol in the directions where the molecular flexibility would be expected to soften the crystal. The other five molecules are all fairly rigid because of the π-electron delocalization in durene, m-dinitrobenzene, resorcinol, and urea and the cage structure of hexamethylenetetramine. While low-energy rotation and torsion of the methyl, amine, and nitro groups will contribute to the computational errors, this appears to be a much smaller effect. Uncertainties in the intramolecular potentials and modeling of such motions, or uncertainties in the placement of hydrogen atoms in a rigid model, will inevitably produce some errors in the calculated elastic constants. Neglect of Temperature and Zero-Point Motion. The calculations presented in this work use second derivatives of the potential energy surface and neglect zero-point and thermal motion. Hence, they strictly only describe the classical, 0 K properties of the crystal. The anharmonicity of the crystal potential energy surfaces leads to thermal expansion and a softening of the crystal. Temperature effects may be responsible for differences of a few percent between the experimental and lattice energy minimized cell lengths. However, the effect of temperature on the elastic constants is considerably greater. In one of the few studies of the temperature dependence of the elastic properties of molecular crystals, the stiffness constants of naphthalene have been measured at several points between 100 and 300 K.105 These values were extrapolated to absolute zero106 to give C11 ) 13.4 GPa, C22 ) 17.3 GPa, and C33 ) 20.8 GPa, which, compared to the room-temperature values of C11 ) 7.8 GPa, C22 ) 10.2 GPa, and C33 ) 11.9 GPa,107 give an estimated mean increase in the stiffness constants of greater than 70% on cooling the naphthalene crystal from room temperature to absolute zero. This contrasts with the minimal increase in the stiffness constants of ionic and metallic crystals. For example, the axial constants of sodium chloride108 and copper15 increase by 18% and 4% from room temperature to 4 K, respectively. If we scale by the melting point to estimate the effects of anharmonicity on the elastic constants between 0 K and room temperature, then durene (mp 352.2 K109) and m-dinitrobenzene (mp 363 K109), which melt at approximately the same temperature as naphthalene (mp 353.5 K109), are expected to soften by approximately the same amount. The hydrogen-bonded crystals β-resorcinol (mp 384 K110), pentaerythritol (mp 542 K109) and urea (mp 408 K109), should soften significantly less than naphthalene, but nevertheless an overestimate of 40-60% would not be unreasonable from comparison with metallic and ionic crystals. Although the elastic constants are generally overestimated, the best rms percent errors in the diagonal terms, about 40% for durene, 40% for m-dinitrobenzene, 40% for resorcinol, 25% for urea, and 50% for hexamethylenetetramine, are not as large as the naphthalene


results would suggest. This may result from the partial incorporation of some of the thermal effects into the empirical exp-6 potential parameters. Although lattice dynamic theory can be used to calculate temperature derivatives of the elastic constants,71,72 ,111-114 it would require an accurate 0 K model intermolecular potential to expect valid comparisons with the experimental elastic constants over a range of temperatures. Although the electrostatic contribution to the potential is readily derived from the molecular wave function, the theory for deriving nonempirical repulsion and dispersion functions from isolated molecule calculations is still in its early development,115,116 despite the problems in representing other minor contributions to the lattice energy. Similarly, the development of methods to incorporate temperature effects into elastic constant calculations requires further experimental measurements for validation over a range of organic systems. The calculation of temperature-dependent elastic constants is likely to prove too expensive for routine use in industrial modeling. Experimental Errors. Although the possible errors due to temperature effects are large, they have to be considered in the context of the errors in the experimental determination. For both urea and pentaerythritol, elastic constants measured at the same (i.e. ambient) temperature differ by more than the combined experimental error estimates. This is indicative of the problems involved with both the experimental measurements and the reduction of raw data (usually either Brillouin scattering or ultrasonic wave velocities) to elastic constants. Some of the problems involved with reducing raw experimental data to elastic constants have been discussed elsewhere.69,117,118 Thus, experimental values may well be in error by more than the reported error estimates. Conclusions The elastic constants of six molecular crystals have been calculated from the curvature of the potential energy surface. An analytic second-derivative method for rigid molecules has been implemented in the crystal structure modeling program DMAREL,59 allowing the use of anisotropic atom-atom potentials. The present work addresses, for the first time, the dependence of the predicted stiffness constants on systematic variations in the model potential. We have found that, for polar molecules, it is essential to represent the electrostatic terms explicitly and preferable to use an anisotropic atom-atom representation. The stiffness of the crystal structure is mainly determined by close contacts; therefore, an accurate representation of the repulsion and dispersion interactions is also important, as shown by treating the methyl and aromatic carbons as different atom types in crystalline durene. The accuracy gained by replacing atomic point charges by an atom-centered multipolar description of the electrostatic interactions is small compared to experimental error. The elastic constants will also be sensitive to the intramolecular potential for flexible molecules. However, when accurate potential models are used, the elastic constants are generally overestimated, which is mostly due to neglect of thermal softening for rigid molecules. This is much more significant for


molecular crystals than for other materials because of the relative weakness of the intermolecular interactions. Nevertheless, the overall agreement with the sparse experimental data shows that these calculations could be used to provide useful predictions of the mechanical properties of organic crystals for the practical design and industrial development of molecular materials. Acknowledgment. We thank Dr. Stuart Walmsley for helpful discussions and the EPSRC for funding of the computing facilities. G.M.D. thanks the University College London Graduate School, the Committee of ViceChancellors and Principals of the Universities of the United Kingdom (Overseas Research Scheme), and the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support. References (1) Payne, R. S.; Roberts, R. J.; Rowe, R. C.; McPartlin, M.; Bashal, A. Int. J. Pharm. 1996, 145, 165. (2) Roberts, R. J.; Payne, R. S.; Rowe, R. C. Eur. J. Pharm. Sci. 2000, 9, 277. (3) Beyer, T.; Day, G. M.; Price, S. L. Manuscript in preparation. (4) Bain, A. M.; El-Korashy, N.; Gilmour, S.; Pethrick, R. A.; Sherwood, J. N. Philos. Mag. B 1992, 66, 293. (5) Gilmour, S.; Pethrick, R. A.; Pugh, D.; Sherwood, J. N. Philos. Mag. A 1993, 67, 855. (6) Manzor, B.; Daly, J. H.; Islam, M. S.; Pethrick, R. A.; Sherwood, J. N. J. Chem. Soc., Faraday Trans. 1997, 93, 3799. (7) Cochran, W. Acta Crystallogr. 1969, A25, 95. (8) Stevens, E. D. Acta Crystallogr. 1974, A30, 184. (9) Harada, J.; Sakata, M. Acta Crystallogr. 1974, A30, 77. (10) Mom, V.; de With, G. Acta Crystallogr. 1978, B34, 2785. (11) Criado, A.; Conde, A.; Marquez, R. Acta Crystallogr. 1985, A41, 491. (12) Day, G. M.; Price, S. L. Elastic Properties of Crystalline Organic Molecules. In The Handbook of Elastic Properties of Solids, Liquids, and Gases; Levy, M., Ed.; Academic Press: New York, 2000; Vol. III. (13) Williams, D. E. J. Chem. Phys. 1966, 45, 3770. (14) Bonadeo, H.; D’Alessio, E. The Refinement of Intermolecular Potentials Using Statical and Dynamical Properties of Molecular Crystals. Application to Chlorinated Benzenes. In Lattice Dynamics and Intermolecular Forces; Corso, L. V., Ed.; Academic Press: New York, 1975. (15) Nelson, D. F. Low-Frequency Properties of Dielectric Crystals, Second and Higher Order Elastic Constants. In LandoltBo¨ rnstein; Numerical Data and Functional Relationships in Science and Technology, New Series. Group III: Crystal and Solid State Physics; Madelung, O., Ed.; Springer-Verlag: Berlin, 1992; Vol. 29, Subvolume A. (16) Allen, F. H.; Kennard, O.; Watson, D. G.; Brammer, L.; Orpen, A. G. J. Chem. Soc., Perkin Trans. 2 1987, S1. (17) Stam, C. H. Acta Crystallogr. 1972, B28, 2630. (18) Sanquer, M.; Ecolivet, C. Elastic Constants in Molecular Crystals: Experiments and Intermolecular Potentials. In Dynamics of Molecular Crystals; Lascombe, J., Ed.; Elsevier: Amsterdam, 1987; p 447. (19) Trotter, J.; Williston, C. S. Acta Crystallogr. 1966, 21, 285. (20) Mitsui, T.; Iio, K. Jpn. J. Appl. Phys. 1980, 19, 2511. (21) Sarkar, S.; Tolapatra, S. K. Fizika (Yugoslavia) 1988, 20, 291. (22) Bacon, G. E.; Lisher, E. J. Acta Crystallogr. 1980, B36, 1908. (23) Koptsik, V. A. Kristallografiia 1959, 4, 219. (24) Semmingsen, D. Acta Chem. Scand. 1988, A42, 279. (25) Srivastava, R. C.; Chakraborty, S. C. J. Phys. Soc. Jpn. 1962, 17, 1767. (26) Srivastava, R. C. Acta Crystallogr. 1962, 15, 1306. (27) Nomura, H.; Higuchi, K.; Kato, S.; Miyahara, Y. Jpn. J. Appl. Phys. 1972, 11, 304. (28) Matsuura, H.; Miyazawa, T. Bull. Chem. Soc. Jpn. 1974, 47, 1143. (29) Swaminathan, B. M.; Craven, B. M.; McMullan, R. K. Acta Crystallogr. 1984, B40, 300.

26


(30) Fischer, G.; Zarembowitch, A. C. R. Seances Acad. Sci., Ser. B 1970, 270, 852. (31) Yoshihara, A.; Bernstein, E. R. J. Chem. Phys. 1970, 77, 5319. (32) Kretchetov, P. M.; Svetlor, V. D.; Teslenko, V. F.; Kitaigorodskii, A. I. Primen. Ul’traaskust. Issled. Veshchestva 1971, 25, 249. (33) Kampermann, S. P.; Sabine, T. M.; Craven, B. M.; McMullan, R. K. Acta Crystallogr. 1995, A51, 489. (34) Haussu¨hl, S. Acta Crystallogr. 1958, 11, 58. (35) Ramachandran, G. N.; Wooster, W. A. Acta Crystallogr. 1951, 4, 431. (36) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular Forces. Their Origin and Determination; Clarendon Press: Oxford, U.K., 1981. (37) Stone, A. J. The Theory of Intermolecular Interactions; Oxford University Press: Oxford, U.K., 1996. (38) Price, S. L. Toward More Accurate Model Intermolecular Potentials for Organic Molecules. In Reviews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; Wiley-VCH: New York, 2000; Vol. 14. (39) Pertsin, A. J.; Kitaigorodskii, A. I. The Atom-Atom Potential Method. Applications to Organic Molecular Solids; SpringerVerlag: Berlin, 1987. (40) Williams, D. E.; Cox, S. R. Acta Crystallogr. 1984, B40, 404. (41) Cox, S. R.; Hsu, L.-Y.; Williams, D. E. Acta Crystallogr. 1981, A37, 293. (42) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J. Phys. Chem. 1996, 100, 7352. (43) Williams, D. E. J. Mol. Struct. 1999, 486, 321. (44) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Chesseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayalla, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. GAUSSIAN 98, Revision A.6; Gaussian Inc.: Pittsburgh, PA, 1998. (45) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (46) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (47) Wiberg, K. B.; Rablen, P. R. J. Comput. Chem. 1993, 14, 1504. (48) Sigfridsson, E.; Ryde, U. J. Comput. Chem. 1998, 19, 377. (49) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (50) Williams, D. E.; Starr, T. L. Comput. Chem. 1977, 1, 173. (51) Gavezzotti, A.; Filippini, G.; Energetic Aspects of Crystal Packing: Experiment and Computer Simulations. In Theoretical Aspects and Computer Modeling of the Molecular Solid State; Gavezzotti, A., Ed.; Wiley: Chichester, U.K., 1997. (52) Rein, R. Adv. Quantum Chem. 1973, 7, 335. (53) Sokalski, W. A.; Poirier, R. A. Chem. Phys. Lett. 1983, 98, 86. (54) Vigne´-Maeder, F.; Claverie, P. J. Chem. Phys. 1988, 88, 4934. (55) Stone, A. J. Chem. Phys. Lett. 1981, 83, 233. (56) Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047. (57) Amos, R. D. with contributions from Alberts, I. L.; Andrews, J. S.; Colwell, S. M.; Handy, N. C.; Jayatilaka, D.; Knowles, P. J.; Kobayashi, R.; Koga, N.; Laidig, K. E.; Maslen, P. E.; Murray, C. W.; Rice, J. E.; Sanz, J.; Simandiras, E. D.; Stone, A. J.; Su, M.-D.; CADPAC6.0: The Cambridge Analytical Derivatives Package: Cambridge, U.K., 1995. (58) Popelier, P. L. A.; Stone, A. J. Mol. Phys. 1994, 82, 411. (59) Willock, D. J.; Price, S. L.; Leslie, M.; Catlow, C. R. A. J. Comput. Chem. 1995, 16, 628. (60) Stone, A. J.; Dullweber, A.; Hodges, M. P.; Popelier, P. L. A.; Wales, D. J. ORIENT, 3.2 ed.; University of Cambridge, Cambridge, U.K., 1996.

Day et al. (61) Buckingham, A. D.; Fowler, P. W. J. Chem. Phys. 1983, 79, 6426. (62) Buckingham, A. D.; Fowler, P. W. Can. J. Chem. 1985, 63, 1985. (63) Buckingham, A. D.; Fowler, P. W.; Hutson, J. M. Chem. Rev. 1988, 88, 963. (64) Price, S. L. J. Chem. Soc., Faraday Trans. 1996, 92, 2997. (65) Ewald, P. Ann. Phys. 1921, 64, 253. (66) Filippini, G.; Gavezzotti, A. Acta Crystallogr. 1993, B49, 868. (67) Gavezzotti, A.; Filippini, G. J. Phys. Chem. 1994, 98, 4831. (68) Drabble, J. R.; Husain, A. H. M. J. Phys. C: Solid State Phys. 1980, 13, 1377. (69) Brose, K.-H.; Eckhardt, C. J. Chem. Phys. Lett. 1986, 125, 235. (70) He, H.; Welberry, T. R. J. Phys. Chem. Solids 1988, 49, 421. (71) Michalski, D.; Swanson, D. R.; Eckhardt, C. J. J. Phys. Chem. 1996, 100, 9506. (72) Michalski, D.; Eckhardt, C. J. J. Phys. Chem. 1997, B101, 9690. (73) Pavlides, P.; Pugh, D.; Roberts, K. J. Acta Crystallogr. 1991, A47, 846. (74) Walmsley, S. H. J. Chem. Phys. 1968, 48, 1438. (75) Walmsley, S. H. Basic Theory of the Lattice Dynamics of Molecular Crystals. In Lattice Dynamics and Intermolecular Forces; Corso, L. V., Ed.; Academic Press: New York, 1975. (76) Sharp, N. D.; Walmsley, S. H. Chem. Phys. Lett. 1994, 222, 546. (77) Catti, M. Acta Crystallogr. 1985, A41, 494. (78) Pavlides, P.; Pugh, D.; Roberts, K. J. Mol. Phys. 1991, 72, 121. (79) Wallace, D. C. Thermodynamics of Crystals; Wiley: New York, 1972. (80) Kulver, R.; Brose, K.-H.; Eckhardt, C. J. Chem. Phys. 1987, 115, 239. (81) Walmsley, S. H. The Calculation of Properties of Molecular Crystals using a Pair Potential Function. In Phonons in Molecular Crystals; Zahlan, A. B., Ed.; University Press: Cambridge, U.K., 1968. (82) Pawley, G. S. Phys. Stat. Sol. 1967, 20, 347. (83) Cerius2; Molecular Simulations Inc.: San Diego, CA, 1997. (84) Sabbah, R.; Tabet, D.; Belaadi, S. T. Thermochim. Acta 1994, 247, 193. (85) Williams, D. E. J. Chem. Phys. 1967, 47, 4680. (86) Nitta, I.; Seki, S.; Momotani, M. J. Chem. Soc. Jpn. Pure Chem. Sect. 1950, 71, 450. (87) Sabbah, R.; Buluku, E. N. L. E. Can. J. Chem. 1991, 69, 481. (88) Barone, G.; Della Gatta, G.; Ferro, D.; Piacente, V. J. Chem. Soc., Faraday Trans. 1990, 86, 75. (89) Ramamoorthy, P.; Krishnamurthy, N. Cryst. Res. Technol. 1997, 32, 525. (90) Ramamoorthy, P.; Krishnamurthy, N. Spectrochim. Acta 1997, A53, 655. (91) McLaghlan, R. D.; Carter, V. B. Spectrochim. Acta 1971, A27, 853. (92) Marzocchi, M. P.; Castellucci, E. J. Mol. Struct. 1971, 9, 129. (93) George, A. R.; Harris, K. D. M.; Rohl, A. L.; Gay, D. H. J. Mater. Chem. 1995, 5, 133. (94) Parlinski, K.; Chapuis, G. J. Chem. Phys. 1999, 110, 6406. (95) Lifson, S.; Hagler, A. T.; Dauber, P. J. J. Am. Chem. Soc. 1979, 101, 5111. (96) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (97) Suponitsky, K. Y.; Tsirelson, V. G.; Feil, D. Acta Crystallogr. 1999, A55, 821. (98) Rousseau, B.; Van Alsenoy, C.; Keuleers, R.; Desseyn, H. O. J. Phys. Chem. 1998, A102, 6540. (99) Keuleers, A.; Desseyn, H. O.; Rousseau, B.; Van Alsenoy, C. J. Phys. Chem. 1999, A103, 4621. (100) Miao, M. S.; Van Doren, V. E.; Keuleers, R.; Desseyn, H. O.; Van Alsenoy, C.; Martins, J. L. Chem. Phys. Lett. 2000, 316, 297. (101) Buergi, H. Personal communication, 2000. (102) DeWit, H. G. M.; Van Miltenburg, J. C.; DeKruif, C. G. J. Chem. Thermodyn. 1983, 15, 651. (103) Belsky, V. K.; Zorkaya, O. N.; Zorky, P. M. Acta Crystallogr. 1995, A51, 473. (104) Mansson, M.; Rapport, N.; Westrum, E. F., Jr. J. Am. Chem. Soc. 1970, 92, 7296.

Elastic Constants for Molecular Organic Crystals (105) Afanasieva, G. K.; Myasnikova, R. M. Kristallografiia 1970, 15, 189. (106) Mirskaya, K. V.; Kozlova, I. E.; Bereznitskaya, V. F. Phys. Stat. Sol. B 1974, 62, 291. (107) Alexandrov, K. S.; Belikova, G. S.; Ryzhenkov, A. P.; Teslenko, V. R.; Kitaigorodskii, A. I. Kristallografiia 1963, 8, 221. (108) Lewis, J. T.; Lehoczy, A.; Bristow, C. V. Phys. Rev. 1967, 161, 877. (109) Beilsteins Handbuch der Organischen Chemie; SpringerVerlag: New York, 1990. (110) Buckingham, J.; Donaghy, S. M. Dictionary of Organic Compounds, 5th ed.; Chapman and Hall: New York, 1982. (111) Della Valle, R. G.; Venuti, E.; Brillante, A. Chem. Phys. 1995, 198, 79.

Crystal Growth & Design, Vol. 1, No. 1, 2001 27 (112) Della Valle, R. G.; Venuti, E.; Brilliante, A. Chem. Phys. 1996, 202, 231. (113) Shpakov, V. P.; Tse, J. S.; Tulk, C. A.; Kvamme, B.; Belosludov, V. R. Chem. Phys. Lett. 1998, 282, 107. (114) Tse, J. S.; Shpakov, V. P.; Belosludov, V. R. J. Chem. Phys. 1999, 111, 11111. (115) Nobeli, I.; Price, S. L. J. Phys. Chem. 1999, A103, 6448. (116) Tsui, H. H. Y.; Price, S. L. CrystEngComm 1999, 7. (117) Huntington, H. B.; Gangoli, S. G.; Mills, J. L. J. Chem. Phys. 1969, 50, 3844. (118) Ahmad, S. F.; Kiefte, H.; Clouter, M. J.; Whitmore, M. D. Phys. Rev. B 1982, 26, 4239.

CG0055070

Elastic Constant Calculations for Molecular Organic Crystals

Recommend Documents