7352
J. Phys. Chem. 1996, 100, 7352-7360
Role of Electrostatic Interactions in Determining the Crystal Structures of Polar Organic Molecules. A Distributed Multipole Study David S. Coombes, Sarah L. Price,* David J. Willock,† and Maurice Leslie‡ Christopher Ingold Laboratories, Chemistry Department, UniVersity College London, 20 Gordon Street, London WC1H 0AJ, England, and EPSRC Daresbury Laboratory, Warrington WA4 4AD, England ReceiVed: February 1, 1996X
The effect of using a realistic model for the electrostatic forces on the calculated structures of molecular crystals is explored by using atomic multipoles derived from an SCF 6-31G** wave function. This was tested on a wide ranging database of 40 rigid organic molecules containing C, H, N, and O atoms, including nucleic acid bases, nonlinear optic materials, azabenzenes, nitrobenzenes, and simpler molecules. The distributed multipole electrostatic model, plus an empirical repulsion-dispersion potential, was able to successfully reproduce the lattice vectors and available heats of sublimation of the experimental room temperature structure in almost all cases. Scaling the electrostatic energy to allow for the effect of electron correlation on the molecular charge density generally improved the lattice energies and the calculated structures to a lesser extent. However, omitting the anisotropic multipole moments usually gave very poor, sometimes qualitatively wrong structures, emphasizing the sensitivity of these crystal structures to the electrostatic forces. We also investigated the sensitivity of the structures to the empirical repulsion-dispersion potential parameters by attempting to optimize these. Since the experimental structures are mainly reproduced to within the errors that could be attributed to the use of static minimization and rigid molecules, it appears that going beyond the atomic charge model to a realistic electrostatic model is a key development in the modeling of the crystal structures of polar and hydrogen-bonded molecules.
1. Introduction The crystal packing adopted by organic molecules reflects the intermolecular forces acting between these molecules, and also determines many physical properties of the solid, ranging from solubility to nonlinear optical activity. Thus, a quantitative understanding of the factors which determine the crystal structure is fundamental to the physical chemistry of condensed phases, as well as being of practical importance in designing new molecular materials with desirable physical properties. It has long been established, by researchers such as Kitaigorodsky,1 that molecules adopt closely packed crystal structures and optimizing the close packing could be sufficient to determine the crystal structure of exceptionally irregularly shaped molecules. However, most molecules can be reasonably close packed in a wide range of structures, requiring an examination of the intermolecular forces to determine the relative stability. The structures of the saturated hydrocarbons, and corresponding lattice energies, can be predicted quite reliably using the atom-atom 6-exp potential to describe the intermolecular forces. This describes the intermolecular potential, U, as the sum of isotropic atom-atom repulsion and dispersion terms, i.e.
U ) ∑ikUik ) ∑ikAικ exp(-BικRik) - Cικ/R6ik
(1)
where atoms i and k are of types ι and κ (C or H), respectively. However, this scheme has not yet been successfully extended to the wide range of organic molecules which includes polar * To whom correspondence should be addressed at the University College London. † Current address: The Leverhulme Centre for Innovative Catalysis, Department of Chemistry, The University of Liverpool, P.O. Box 147, Liverpool L69 3BX, England. ‡ EPSRC Daresbury Laboratory. X Abstract published in AdVance ACS Abstracts, April 15, 1996.
S0022-3654(96)00333-4 CCC: $12.00
and hydrogen-bonding functional groups.1 Recently, Filippini and Gavezzotti have determined a set of 6-exp parameters for modeling crystal structures containing nonpolar2 and monofunctional hydrogen-bonding compounds.3 These potentials are remarkably successful at predicting the lattice energies of the molecular crystals, but the inability to accurately reproduce the experimental crystal structures was partially ascribed to the neglect of electrostatic anisotropy. The atom-atom potential can be extended to represent the electrostatic forces by including the interactions between point charges on the atomic centers. Such an isotropic atom-atom 6-exp-1 model is widely used in molecular modeling, for its simplicity. The results of crystal modeling studies with such models generally reinforce the conclusion of Desiraju4 that this simple model of intermolecular interactions “gives an approximate description but does not stress adequately many structure-defining interactions”. It has already been found that the errors in the atomic charge model have a significant effect on the calculated crystal structures of the azabenzenes5 by Williams and Weller. They found it necessary to introduce additional lone pair charge sites off the nitrogen atoms, which also significantly improved the fit to the electrostatic potential around the molecules, in order to give a satisfactory minimum energy crystal structure. The packing symmetry which gives the lowest lattice energy for a diatomic molecule X2 is also sensitive to the form of the electrostatic model, and can depend on whether the same total quadrupole moment is represented by atomic dipoles or atomic quadrupoles.6 The electrostatic forces between molecules can be modeled far more accurately by representing the molecular charge distribution by sets of point multipoles (charge, dipole, quadrupole, etc.) on every atomic site. Such distributed multipole electrostatic models automatically include the effects of lone pair and π electron density on the intermolecular forces. They © 1996 American Chemical Society
Crystal Structures of Polar Organic Molecules
J. Phys. Chem., Vol. 100, No. 18, 1996 7353
Figure 1. Molecules whose crystal structures were studied. Along with the common name of the molecule, the crystal structure is denoted by the refcode of the entry in the Cambridge Structural Database.15
are being increasingly widely used to model gas phase interactions, and thereby understand the structures found in van der Waals complexes of polar7 and aromatic molecules,8 and for predicting the structures of protein side chains9 and nucleic acid base pairs.10 This appears to be successful since the electrostatic
forces generally dominate the anisotropy of the intermolecular interactions.11 Since the use of distributed multipole, or other accurate, electrostatic models has provided a major qualitative improvement in our ability to model van der Waals complexes, we expect the use of such models to likewise improve our ability
7354 J. Phys. Chem., Vol. 100, No. 18, 1996
Coombes et al.
TABLE 1: Calculated Molecular Crystal Structures, Using Distributed Multipole Based Model Potentials crystal refcode
errors in calculated structures (calcd - exptl) [first line EST+0.9DMA, second line FIT+0.9DMA for molecules containing polar hydrogen atoms]
experimental structure a/Å
b/Å
c/Å
β/dega
∆a/Å
∆b/Å
∆c/Å
∆β/dega
∆θ/deg
∆x/Å
Fb
URICAC
14.464
7.403
6.208
65.10
GEBTUC
7.102
9.759
10.387
58.85, 67.64, 72.00
0.436 0.284 -0.294
0.029 0.127 0.241
0.030 -0.074 0.502
2.849 1.467 3.786
0.318 0.228 0.356
23.50 16.69 78.71
-0.301
0.123
0.486
4.100
0.362
74.65
0.226 0.164 0.293 0.318 0.016 0.011 -0.069 -0.277 0.269 0.209 -0.214 -0.096 -0.336 -0.429 0.181 0.084 -0.331 -0.278 0.002
0.2803 0.1933 -0.534 -0.752 0.006 -0.129 -0.003 -0.073 0.199 0.138 0.186 0.059 0.451 0.451 -0.180 -0.122 0.534 0.549 -0.062
0.0439 0.0497 -0.066 -0.086 0.030 0.009 0.042 0.068 0.238 0.084 0.000 0.002 0.061 0.100 0.191 0.064 0.955 0.696 -0.062
-1.36 -1.64 1.35, -3.09, -2.13 1.70, -2.99, -2.02 3.16 3.19
2.932 1.827 4.592 4.684 2.142 2.248 2.569 3.832 2.185 2.107 7.533 6.963 3.363 3.681 2.381 3.118 11.099 11.545 0.615
0.275 0.247 0.341 0.342 0.083 0.087 0.216 0.247 0.229 0.196 0.168 0.138 0.335 0.365 0.264 0.240 0.583 0.478 0.380
29.85 23.30 63.93 80.88 6.75 7.07 7.83 18.10 39.15 22.21 22.82 14.81 43.77 58.27 43.79 25.83 383.20 276.62 17.17
-0.010
-0.135
-0.135
0.631
0.467
34.02
-0.031
0.435
-0.280
4.842
0.185
67.72
-0.134
0.474
-0.334
6.555
0.214
98.36
-0.362 -0.443 0.34 -0.110 0.254 0.480 0.351 0.132
-0.141 -0.231 0.43 0.202 -0.074 -0.241 -0.282 -0.100
-0.073 -0.170 -0.40 -0.154 -0.184 -0.074 -0.134 -0.029
1.829 2.079 26.704 2.254 2.013 3.140 2.709 4.741
0.127 0.157 0.324 0.141 0.114 0.266 0.254 0.189
20.49 37.48 398.72 25.03 23.18 47.64 49.33 32.52
0.057
-0.121
-0.067
2.877
0.156
20.76
0.060 0.011 -0.1304 -0.1500 -0.20 -0.27 -0.075 -0.167 0.57 0.47 0.285 0.491 0.096 0.060 -0.204 -0.317 -0.709 -0.684
0.065 0.056 -0.0066 -0.0187 -0.09 -0.19 0.091 0.033 -0.35 -0.33 -0.389 -0.518 -0.248 -0.291 -0.031 -0.031 0.319 0.309
2.111 2.044 0.598 0.577 15.144 13.783 0.000 0.000 3.683 2.993 4.274 6.272 2.420 2.619 4.439 4.338 5.766 5.587
0.202 0.179 0.109 0.135 0.404 0.382 0.046 0.093 0.399 0.344 0.385 0.532 0.092 0.109 0.271 0.312 0.450 0.442
14.14 11.26 4.01 5.98 265.32 218.41 7.50 18.95 73.40 56.03 96.49 175.06 13.32 15.91 24.97 29.76 130.35 119.86
URACIL
11.938
12.3760
3.6552
AZURAC10
4.875
17.611
5.022
ALOPUR
3.683
14.685
10.318
CYTSIN01
13.044
9.496
3.814
CYURAC10
7.900
6.732
11.951
15.553
9.374
3.664
PYRZOL02
8.232
12.840
7.054
IMAZOL13
7.732
5.458
9.779
117.26
INDAZL
7.570
5.760
7.700
119.5
FITXIP
6.257
6.129
6.129
PURINE
DAMTRZ20
10.658
10.837
4.339
120.90
97.47
130.67
79.508, 90.000, 90.000
90.00, 90.00, 118.83
MELAMI04
10.606
7.495
7.295
112.26
TETRAZ PRMDIN PYRAZI JOWWIB
5.23 11.698 9.316 14.199
5.79 9.493 3.815 4.844
6.63 3.806 5.911 14.258
115.25
TRPHAM
5.027
5.355
7.165
BZAMID01
5.607
5.046
16.557
10.2453
FORMAM
3.69
UREAXX09
5.645
BICWEP
ONITAN
15.45
22.053
105.98 103.02, 100.90, 92.17
90.66
5.9400
91.840
9.18
6.87
98.0
5.645
4.704
10.01
8.57
107.4
MNIANL01
6.499
19.369
5.084
NANILI21
12.337
6.037
8.597
91.42
CUZDEG
3.782
27.729
7.152
95.41
TNIOAN
6.137
9.217
15.323
99.67
-0.104 -0.111 -0.082 -0.168 0.12 0.08 -0.075 -0.167 -0.35 -0.34 0.264 0.337 0.045 0.086 0.055 0.057 0.323 0.287
-2.15 -2.04
2.90 2.36
3.90 3.32 7.7 5.9 -0.787, 0.000, 0.000 -1.600, 0.000, 0.000 0.00, 0.00, 0.84 0.00, 0.00, 1.70 1.38 1.26 8.66 1.30 0.39 2.56, 2.51, 0.04 2.11, 1.67, -1.32 -1.99 -1.74 -0.923 -0.893 13.2 11.6 0.2 -0.2
0.14 -0.24 3.15 3.39 3.24 3.38
notesc
FS
A∆T FS ∆T FSL A∆T FSL FS
N sp3? FS
N sp3?
N sp3? FS FS
A∆T
N sp3? FSL N sp3? FSL N sp3? FSL FS N sp3? FSL
Crystal Structures of Polar Organic Molecules
J. Phys. Chem., Vol. 100, No. 18, 1996 7355
TABLE 1 (Continued) crystal refcode TATNBZ
errors in calculated structures (calcd - exptl) [first line EST+0.9DMA, second line FIT+0.9DMA for molecules containing polar hydrogen atoms]
experimental structure a/Å
b/Å
c/Å
β/dega
∆a/Å
∆b/Å
∆c/Å
∆β/dega
∆θ/deg
∆x/Å
Fb
9.010
9.028
6.812
108.58, 91.82, 119.97
-0.057
-0.077
-0.584
12.835
1.171
521.58
-0.107
-0.123
-0.597
12.964
1.154
526.35
-0.309 -0.415 -0.27 -0.003 -0.152 -0.333 0.1064 -0.057 0.051 0.040 0.165 0.40
0.275 0.341 0.16 -0.040 -0.143 -0.117 -0.0669 0.006 -0.101 -0.109 -0.099 -1.93
-0.175 -0.290 -0.06 -0.359 -0.021 0.164 0.1579 -0.117 -0.026 -0.062 -0.224 1.10
-15.16, -6.23, 0.09 -15.27, -6.17, 0.08 0.39 0.27
6.238 7.006 1.560 3.933 3.945 4.771 3.692 6.843 1.934 1.637 1.957 25.290
0.368 0.444 0.186 0.158 0.157 0.135 0.121 0.083 0.058 0.068 0.080 0.799
51.68 81.35 12.04 34.01 9.00 29.43 14.29 24.30 6.39 7.66 23.50 897.55
BAZGOY
21.822
5.867
8.386
TNBENZ10 ZZZFYW01 DNBENZ10 DNITBZ11 OBNZQU BNZQUI NAPHTD ANTCEN14 PHENAN08 BENZEN02
9.78 7.945 13.257 11.137 6.3137 7.055 8.266 8.553 8.441 7.46
26.94 12.975 14.048 5.461 5.7985 6.795 5.968 6.016 6.140 9.66
12.82 7.421 3.806 5.684 6.8865 5.767 8.669 11.172 9.438 7.03
101.01 111.88 92.22 101.47 122.92 124.596 97.96
-2.05 0.24 -2.67 -1.33 -1.640 -3.14
notesc
N sp3?
∆T
∆T
a
For monoclinic systems, β and ∆β are given. No entry implies a tetragonal or orthorhombic crystal system. For other systems, cell angles are given in the order R, β, γ and cell angle errors in the order ∆R, ∆β, ∆γ. b F, the structural drift factor, is defined2 as F ) (∆θ/2)2 + (10∆x)2 + (100∆a/a)2 + (100∆b/b)2 + (100∆c/c)2 + ∆R2 + ∆β2 + ∆γ2, where rms values for ∆θ (the net rotation of the molecule (deg)) and ∆x (the net translation of the center of mass of the molecule (Å)) are used when there is more than one independent molecule in the asymmetric unit cell.c N sp3? denotes molecules in which the NH2 group is probably significantly sp3 hybridized as the C-N bond to the aromatic ring is longer than 1.34 Å. Thus, there is additional uncertainty in the hydrogen atom position. ∆T indicates that there is a low-temperature crystal structure available which suggests that the molecule has a particularly large thermal expansion coefficient (>0.03 Å3 K-1, >0.05 Å3 K-1 for hydrocarbons). A∆T indicates that a cell edge expands on cooling. This suggests that the error in using static minimization will be particularly large for these molecules. However, this assessment can only be made for a minority (12) of these molecules. FS denotes structures used in the fitting of the reported Hp parameters, with FSL denoting that the lattice energy was also used.
to model crystal structures, which also contain molecules in van der Waals contact. Indeed, atomic multipole moments derived from X-ray data have already been used to account for the crystal packing characteristics of amides12 and carboxylic acids.13 Thus, we wish to test whether the introduction of realistic, distributed multipole, electrostatic models into the computer simulation of molecular crystal structures produces a qualitative improvement in our ability to model a wide range of polar and hydrogen-bonded molecular crystals. In this paper, we use an accurate model for the intermolecular electrostatic interactions, derived by a distributed multipole analysis (DMA)14 of an ab initio molecular charge distribution, to study the crystal packing of polar molecules ranging from nucleic acid bases to nitrobenzenes. The wave function (and hence DMA) automatically represents the variations in the atomic charge distributions with hybridization and chemical environment. Otherwise, we follow the traditional approach to modeling molecular crystal structures, by minimizing the static lattice energy of rigid molecules, and by using empirical 6-exp models to represent the other contributions to the intermolecular potential. 2. Methods The polar organic molecules considered in this study ranged from nucleic acid bases and amides to nitrobenzenes and nitroanilines, to cover examples of some simple organic pharmaceutical, nonlinear optic, and explosive materials. The molecules are shown in Figure 1. They were chosen to sample various combinations of aromatic rings, heterocyclic rings, and amine (NH2), nitro (NO2), amide (C(O)NH) groups, and so some of the molecular crystals include NsH‚‚‚OdC and NsH‚‚‚N hydrogen bonds and π-π interactions. Relatively small molecules which could be reasonably modeled as rigid, and whose shapes were unlikely to determine the crystal packing, were chosen. The choice was made, with some degree of
arbitrariness, from molecules where a reasonable quality, approximately room temperature crystal structure was found in the Cambridge Structural Database,15 and includes quasi-layered structures, but excludes most well-characterized polymorphs. Their structures were taken directly from the experimental crystal structures, with the hydrogen atoms positioned to give a standard bond length16 of 1.08 Å for CsH and 1.01 Å for NsH along the experimental bond direction. The calculations used an atom-atom model intermolecular potential which comprised an ab initio based distributed multipole model for the electrostatic contribution and an empirically fitted 6-exp repulsion dispersion. The electrostatic model is the main novel feature and focus of this work, as atomcentered multipole expansions have only been previously tested for their ability to reproduce the crystal structures of some explosive molecular crystals.17 The DMA electrostatic model was obtained from the SCF wave function of the isolated molecule, calculated for a 6-31G** basis set,18 within the program suite CADPAC.19 This gave sets of atomic multipoles, up to hexadecapole, on every atomic site. Thus, up to 25 independent multipole components, depending on the molecular symmetry, were used to represent the electronic charge distribution around each atom. The electrostatic interactions were calculated using all terms up to R-5 in the multipole expansion,20 using Ewald summation methods for the charge-charge, charge-dipole, and dipoledipole contributions to the lattice energy,21 and direct summation over all atoms in molecules whose centers were within 20 Å for the rest of the electrostatic contribution to the lattice energy. Thus, the main error in the evaluation of the electrostatic contribution to the lattice energies, excluding the effect of penetration into the molecular charge distributions, is the accuracy of the ab initio wave function. It is well known that the neglect of electron correlation inherent in an SCF wave function leads to an overestimate of the molecular dipole
7356 J. Phys. Chem., Vol. 100, No. 18, 1996
Coombes et al.
moment by 12-15%.22 A comparison of experimental dipole moments with those calculated from a 6-31G** wave function for eight small molecules, by Cox and Williams,23 led to the scaling factor of 0.9 which is widely used as a “fudge factor” in molecular modeling. This factor is also justified by a study of the effect of basis set and electron correlation on the electrostatic interactions of a dipeptide calculated from DMAs.24 These results also suggest that the SCF overestimate of the electrostatic interaction is predominantly a scaling effect for reasonable quality basis sets, so the anisotropy of the electrostatic potential around the molecule should be realistic. Thus, three electrostatic models were tested, the full DMA of the SCF 6-31G** wave function (DMA), its isotropic component (i.e., only the charge terms) to investigate the importance of the anisotropic multipoles (Q), and the full set of multipoles multiplied by the scaling factor of 0.9 (0.9DMA). The last model should be the best approximation to the actual electrostatic forces between the molecules. Ideally, the other terms in the intermolecular potential should also be derived from the molecular wave functions. However, although methods of deriving repulsion and dispersion potentials from the molecular wave function are being developed, unlike the electrostatic contribution, they cannot yet be routinely applied to organic molecules.25 Hence, the sum of all other contributions to the intermolecular potential was assumed to have the isotropic atom-atom form (1), and depend only on the separation Rik of the atoms in different molecules, and the types of the atoms concerned. The types considered were carbon (C), nitrogen (N), oxygen (O), hydrogen bonded to carbon (Hc), and hydrogen bonded to nitrogen (Hp), to keep the number of types to a minimum. The distinction between polar hydrogen atoms Hp and nonpolar Hc is necessary to model the smaller effective van der Waals radius of protons involved in hydrogen bonding. This reflects the major difference in the electron density associated with the two types of hydrogen atoms. We are effectively assuming that the variations in the charge distributions associated with different hybridization and bonding environments of nitrogen, oxygen, and carbon atoms need only be represented in the electrostatic model, and that all other contributions to the intermolecular potential are represented by the empirical transferable repulsion-dispersion potential. The empirical potential repulsion-dispersion parameters were taken from Williams and Cox26 for C, Hc, and N, derived from fitting to azahydrocarbons, with the O parameters from the fitting to oxohydrocarbons.27 These potentials were fitted empirically, in conjunction with potential derived charges, including lone-pair charge sites on the azabenzenes. These potentials were not designed for hydrogen-bonded crystals, and so the starting point for the polar hydrogen potentials was taken from the O‚‚‚Hp potential fitted to intermolecular perturbation theory calculations of the exchange-repulsion, penetration and dispersion interaction between formamide and formaldehyde in the NsH‚‚‚OdC hydrogen-bonding region.28 The heteroatomic parameters were fixed using the traditional combining rules
Aικ ) (AιιAκκ)1/2, Bικ ) (1/2)(Bιι + Bκκ), Cικ ) (CιιCκκ)1/2 (2) The combining rule for the Cικ dispersion coefficients has some physical justification at the level of the Unso¨ld approximation to the R-6 dispersion coefficient, but the other rules, though widely used, are poorly justified and known to be limited in their accuracy.29 These repulsion-dispersion potentials have not, to our knowledge, been tested for molecules containing both N and O, or including hydrogen bonding, let alone mixed
functional groups. This composite repulsion-dispersion potential will be denoted EST to emphasize that it has been put together from among the most appropriate literature sources but is nevertheless being used outside the realm for which the individual parameters were derived and tested. Attempts to improve on this estimate by fitting are described later. A number of the structures investigated contain exocyclic NH2 groups, where the C-N bond length is longer than 1.34 Å, suggesting a significant degree of sp3 hybridization and thus further uncertainty in the direction of the bond to the hydrogen atoms. Thus, we also investigate the effect of altering the hydrogen atom positions by making the NH2 groups approximately coplanar with the ring systems of the molecules concerned. For benzene and indazole, we also investigate the effect of crudely representing the shift of the hydrogen electron density into the bond30 by placing their interaction sites at 1.01 Å along the C-H bonds, i.e., a distance of 0.07 Å from the hydrogen nuclear positions used in determining the wave function. The anisotropy of the atom-atom electrostatic interactions required a purpose-written crystal structure analysis code DMAREL.31 The code calculates the noncentral forces and torques arising from all the terms in the multipole expansion, and calculates the net force, the torque about the center of mass of each molecule in the unit cell, the strains on the lattice, and the second derivatives of the lattice energy. A modified Newton-Raphson procedure is used to minimize the lattice energy with respect to the six elements of the strain matrix which describes the changes in the shape of the unit cell, and the six Cartesian components of the center of mass vector and orientational vector of each molecule in the primitive unit cell. Since Cartesian coordinates are used, the symmetry of the experimental crystal structure, used as a starting point, is not a constraint on the minimization. 3. Results 3.1. Variations in the Electrostatic Model. The errors in the calculated crystal structures for the theoretically best justified model EST+0.9DMA are given in Table 1. Half this diverse range of structures are reproduced with rms errors in the lattice parameters of less than 2%, including examples from each type of functional group. A high proportion of those structures calculated with errors of 2-5% include exocyclic NH2 groups where the C-N bond lengths suggest a significant degree of sp3 hybridization, and thus an uncertainty in the hydrogen atom positions, which limits the accuracy of the calculated structure (see section 3.4 and Table 4). Only three structures are calculated to have a qualitatively wrong structure, and these are among the simplest molecules: benzene, s-tetrazine and indazole (see section 3.4). The effect of removing the 0.9 scaling factor, i.e., omitting any consideration of the effect of electron correlation on the wave function, and then the further removal of the anisotropic atomic multipoles is shown in Table 2. The scaling of the electrostatic energy by 0.9 × 0.9 ) 0.81 has a noticeable effect, though on average it affects the lattice vectors by less than 1%. The scaled results are generally better, in accord with the belief that the electrostatic forces are generally being calculated more accurately, but there are many exceptions. In contrast, the omission of the anisotropic multipole moments has a major effect on most structures, with only cyanuric acid and FITXIP giving less than 2% rms error in the lattice parameters. Many structures are very poorly reproduced with several diverse molecules, such as m-dinitrobenzene, allopurinol, formamide, uric acid, uracil, and o- and p-nitroaniline, being calculated to have a qualitatively wrong structure.
Crystal Structures of Polar Organic Molecules
J. Phys. Chem., Vol. 100, No. 18, 1996 7357
TABLE 2: Sensitivity of the Calculated Crystal Structure to the Intermolecular Potential rms % error in the lattice parametersa variations in electrostatic model crystal
EST+ 0.9DMA
EST+ DMA
EST+Q
URICAC GEBTUC URACIL AZURAC10 ALOPUR CYTSIN01 CYURAC10 PURINE PYRZOL02 IMAZOL13 INDAZL FITXIP DAMTRZ20 MELAMI04 TETRAZ PRMDIN PYRAZI JOWWIB TRPHAM BZAMID01 BICWEP FORMAM UREAXX09 ONITAN MNIANL01 NANILI21 CUZDEG TNIOAN TATNBZ BAZGOY TNBENZ10 ZZZFYW01 DNBENZ10 DNITBZ11 OBNZQU BNZQUI NAPHTD ANTCEN14 PHENAN08 BENZEN02
1.8 3.9 1.8 4.0 0.3 0.7 2.8 1.4 3.1 2.6 9.3 0.7 4.4 2.3 6.7 2.7 2.6 3.5 1.9 1.3 0.8 2.3 1.7 4.2 5.1 1.9 1.0 5.5 5.0 3.1 1.6 2.8 0.9 2.7 1.8 1.3 1.1 1.1 2.0 15.0
1.3 3.3 1.4 4.2 0.6 2.2 1.7 0.5 3.3 2.5 5.6 1.5 4.4 3.2 7.9 2.5 2.8 3.1 1.5 0.9 1.1 3.2 1.6 3.6 6.1 2.1 1.4 5.5 5.3 4.5 1.6 2.9 1.2 2.4 1.5 1.5 1.8 2.0 2.6 14.6
20.8 5.4 25.4 8.0 11.6 4.4 1.0 5.9 0.9(nc) 6.7 6.1 1.7 5.8 2.5 24.5 5.5 9.5 8.3 6.0 4.4 4.1 12.9 3.0 13.2 5.4 27.2 38.0(nc) 4.3 9.0 2.5 3.3 2.8 40.3 3.8 3.3 8.1 3.2 3.9 5.2 7.9
other variations in potential FIT+ 0.9DMA
publishedb
1.7 3.7 1.4 4.6 0.5 1.7 2.0 0.5 3.7 1.5 7.9 1.6 5.2 3.3 7.4 LP5 3.0 LP5 2.6 LP5 3.7 1.6 1.2 1.0 2.6 2.2 3.7 6.8 2.1 1.1 5.2 5.2 4.1
4.3 LJ35 1.2 LJ35
2.0 EH36
4.127 1.326 0.926 1.226 1.726
The rms % error is calculated as 3δ2 ) (100∆a/a)2 + (100∆b/b)2 + (100∆c/c)2 except for the higher symmetry tetragonal cell of UREAXX09 and the pseudomonoclinic cell of FITXIP, where the average is taken over the two independent cell edges. b Values of δ from other published work, where available, based on isotropic 6-exp atom-atom potentials, usually with atomic charges. Significant variations in the potential form are denoted: LP, additional lone pair charge sites included; LJ, 12-6 repulsion-dispersion potential; EH, empirical hydrogen bond potential with no electrostatic model. nc denotes minimizations which were not properly converged. a
3.2. Lattice Energies. The lattice energies predicted by the different models differ significantly because the electrostatic contribution to the lattice energy at the experimental structure is very significant for all but the hydrocarbons, as shown in Table 3. The full DMA calculations always have a lower lattice energy, as the electrostatic contribution is always attractive. The 0.9 scaling factor reduces the lattice energy by less than 5 kJ mol-1 for aromatic hydrocarbons to over 20 kJ mol-1 for multiply hydrogen-bonded systems, so the reliability of the electron correlation fudge factor is important in the estimated lattice energies. The charge-only calculations often predict a very different lattice energy from the full DMA. This can be attributed to the effects of the higher multipole moments, and is particularly marked for the molecules containing amide and nitro groups. However, the effect is naturally dependent on
the scheme used to partition the charge density, and so much better lattice energies should be obtained by using potential derived atomic charges. The comparison of the calculated lattice energies with the experimental sublimation energy is problematic, because of the neglect of zero point energy and the high experimental error.32 An analysis of these problems caused Kitaigorodsky to conclude that discrepancies up to 3-4 kcal mol-1 (12-17 kJ mol-1) should not cause any concern when judging the quality of a potential model. Many lattice energies are being predicted within this limit by either DMA model. However, the EST+0.9DMA model provides the best predictions of the lattice energies of the non-hydrogen-bonded molecules, but generally underestimates the lattice energies of the hydrogen-bonded molecules. This is consistent with the scaling of the electrostatic interactions being required to allow for an overestimate of the electrostatic energy by an SCF wave function, but the estimated polar hydrogen parameters do not include the attractive polarization and charge-transfer contributions and thus underestimate the well depth of a hydrogen bond. 3.3. Optimization of the Repulsion-Dispersion Parameters. The repulsion-dispersion parameters were taken from various literature sources, and not developed for many of the interactions being sampled in this diverse range of crystal structures, and so we attempted to optimize them for use in conjunction with a DMA electrostatic model. Various fitting schemes were tried, based either on minimizing the forces, torques, and strains or alternatively on minimizing the first step displacement that would occur on the energy minimization. Various subsets of up to 20 molecular crystal structures and associated lattice energies were used in the fitting. Although most fits sought to optimize the preexponential repulsion parameters (Aιι), several variants, such as relaxing the combining rules, or separating the N atoms into different types, were attempted. The results can be categorized into one of two classes. If the fitted parameters changed significantly from the starting values of EST, then the calculated structures were significantly poorer. (This usually resulted from an attempt at a poorly defined fitting problem.) On the other hand, most wellbehaved fits resulted in parameters that differed little from the starting point and resulted in calculated structures that, while different for most individual molecules, were overall of the same quality as those reported in Table 2 for EST+0.9DMA. The main deficiency in EST+0.9DMA is that the calculated polar hydrogen parameters had not been derived in conjunction with the other empirical parameters, and did not include the polarization and charge-transfer contributions to the hydrogenbonding energy. Optimization of the Hp parameters to the dataset of thirteen hydrogen-bonded structures and six lattice energies (defined in Table 1) produced a decrease in AHpHp from 7017.3 to 5029.6 kJ mol-1 and an increase in CHpHp from 16.4 to 21.5 kJ mol-1 Å6. The structures calculated with the new parameters (FIT+0.9DMA) show sensitivity to the Hp repulsiondispersion parameters (Table 1), but the overall quality is the same. The deeper hydrogen-bonding potential well improves the agreement with the lattice energies somewhat, but they are still generally underestimated. The predicted lattice energies (Table 3) are in the range that should not give concern in comparison with the heats of sublimation. However, the uncertainty in the electrostatic scaling factor and the probable inability of this repulsion-dispersion model to completely absorb the charge-transfer and polarization contributions, combined with the likely experimental error in the known sublimation energies, reduce the reliance that can be put on the absolute values of the calculated lattice energies.
7358 J. Phys. Chem., Vol. 100, No. 18, 1996
Coombes et al.
TABLE 3: Lattice Energiesa lattice energies/(kJ mol-1) error (∆Hsub - Ulatt) in calculated energy at relaxed structure
at experimental structure
model
∆Hsub/ (kJ mol-1)
EST+0.9DMA
EST+DMA
EST+Q
FIT+0.9DMA
0.9DMA
EST+0.9DMA
GEBTUC URACIL AZURAC10 CYTSIN01 CYURAC10 PYRZOL02 IMAZOL13 INDAZL MELAMI04 PRMDIN PYRAZI BZAMID01 FORMAM UREAXX09 ONITAN MNIANL01 NANILI21 TNIOAN TATNBZ TNBENZ10 ZZZFYW01 DNBENZ10 DNITBZ11 BNZQUI NAPHTD ANTCEN14 PHENAN08 BENZEN02
-158.1 -131.0 -141.0 -155.0 -133.0 -74.0 -83.0 -91.1 -123.3 -49.9 -56.3 -101.7 -71.7 -98.6 -90.0 -108.3 -107.0 -125.3 -168.2 -107.3 -87.9 -87.0 -96.2 -62.8 -72.5 -97.8 -90.9 -44.4
-24.7 -23.9 -44.5 -23.2 -21.8 -11.0 -14.1 -14.1 7.6 7.3 -2.8 -6.8 -5.4 -20.3 -8.9 -13.9 -11.7 -5.5 -7.0 0.0 12.8 3.1 4.0 11.2 3.1 6.7 10.2 5.0
-1.1 -6.0 -29.3 1.5 -0.5 0.1 -2.5 -6.2 36.0 12.1 1.2 5.1 9.1 -2.8 -3.1 -3.6 -1.6 4.1 11.9 8.4 20.6 8.3 11.8 18.3 5.1 9.6 12.7 7.4
31.6 33.4 21.4 15.3 92.3 -26.8 -13.9 -23.4 -10.7 10.3 -1.6 19.2 38.9 26.7 16.7 -2.3 18.7 105.4 133.0 119.5 41.6 62.6 68.1 58.4 -3.3 -3.0 1.0 -1.8
-18.6 -19.7 -40.2 -15.8 -14.5 -7.8 -11.1 -12.0 21.8
-110.5 -85.6 -70.5 -105.7 -103.3 -46.4 -51.4 -26.3 -97.2 -20.9 -17.0 -51.2 -63.1 -75.1 -23.6 -37.6 -41.7 -40.6 -60.6 -35.2 -32.3 -21.5 -32.3 -30.7 -8.4 -11.8 -10.0 -6.6
-129.7 -102.7 -93.0 -128.5 -105.1 -61.5 -67.2 -73.6 -126.7 -56.0 -52.6 -94.1 -61.8 -77.5 -79.3 -92.6 -94.6 -118.5 -154.0 -105.7 -98.1 -88.5 -97.5 -71.6 -74.7 -102.9 -99.0 -47.7
-2.5 -1.0 -13.8 -7.4 -11.0 -9.0 -3.8 -0.1
a Experimental heats of sublimation (∆H ) taken from ref 37, except PRMDIN.38 Where more than one experimental determination has been sub made, the most recent or largest result was usually used, since experimental heats of sublimation are more easily underestimated.39 The lattice energy was evaluated by direct summation to 20 Å for all terms except the charge-charge, charge-dipole, and dipole-dipole electrostatic contribution, which were evaluated by Ewald type summation.
TABLE 4: Effect of Hydrogen Atom Positions on Calculated Crystal Structures
crystal refcode
F
δ/%
δ/% for X-H direction from experimental structure
154.15
6.3
4.4
-9.9
22.79 75.68 47.33 9.65 12.16 129.39 32.80
2.2 4.2 3.4 1.4 1.8 5.5 2.9
2.3 3.5 4.2 5.1 1.9 5.5 3.1
-17.1 -2.3 0.3 -5.0 -0.6 1.1 -4.7
134.76 384.39
3.7 9.1
15.0 9.3
1.8 2.5
errors in EST+0.9 DMA for altered H atom positions ∆a/Å
∆b/Å
∆c/Å
DAMTRZ20
0.014
0.720
-0.373
MELAMI04 JOWWIB ONITAN MNIANL01 NANILI21 TNIOAN BAZGOY
0.077 0.620 -0.29 0.002 0.041 0.320 0.170
-0.277 -0.282 0.46 0.008 0.092 -0.706 -0.252
0.018 0.070 -0.28 -0.125 -0.234 0.330 0.207
BENZEN02 INDAZL
-0.47 -0.376
-0.12 0.493
0.03 0.942
∆β/deg
∆θ/deg
∆x/Å
Use of Idealized sp2 NH2 Group 0.00, 10.057 0.331 0.00, 0.15 -0.44 5.354 0.109 2.01 5.987 0.306 -0.4 2.789 0.313 3.655 0.052 0.20 2.508 0.082 3.23 5.800 0.449 -0.12 3.777 0.202 Use of Foreshortened C-H Bondb 18.972 0.172 8.2 11.193 0.611
∆Ulatta/ (kJ mol-1)
a ∆Ulatt ) Ulatt(X-H direction from experimental structure) - Ulatt(new H positions). b The interaction site was placed 1.01 Å along the C-H bond, a distance of 0.07 Å from the H nuclear position used in determining the wave function.
3.4. Exceptions and Comparison with Published Work. A detailed, realistic electrostatic model with an isotropic atomatom repulsion-dispersion potential is able to account for the crystal structures of a wide range of organic molecular crystals to within a few percent error in the lattice parameters. The effects of altering the hydrogen atom positions for those structures containing exocyclic NH2 groups which may be partially sp3 hybridized are shown in Table 4. These results show that some structures are very sensitive to the assumed static position of the hydrogen atoms, while others are fairly insensitive. Thus, the unknown dynamic behavior of sp3 amine
groups within crystals limits the accuracy to which the experimental structures can be reproduced. The structures which are poorly reproduced by this potential scheme are among the simplest: benzene, s-tetrazine, and indazole. The benzene crystal structure undergoes a molecular rotation of 25° on minimization with a gain of only 1.7 kJ mol-1 in lattice energy. With the hydrogen interaction sites shifted by 0.07 Å toward the carbon atoms along the C-H bonds, an rms error of 3.7% was obtained (see Table 4), demonstrating the sensitivity of this structure to the anisotropy of the H atom interactions. The calculated crystal structure of benzene will
Crystal Structures of Polar Organic Molecules also be particularly dependent on the simulation method, given the large thermal expansion and dynamic behavior. The packing of the carbon aromatic ring in indazole suggests that it may be poorly calculated for similar reasons, but it was found to be much more sensitive to the scaling of the electrostatic forces than the position of the hydrogen interaction site. The crystal structure of s-tetrazine is also probably unusually sensitive to the anisotropy of the repulsion potential, as it has unusual N‚‚‚N close contacts and has been reproduced accurately only by an anisotropic repulsion potential.33 A comparison of these calculated crystal structures with the limited number of other published ones is given in Table 2. The structures calculated using FIT (or EST)+0.9DMA are very comparable for all structures except benzene and 1,3,5-triamino2,4,6-trinitrobenzene, despite the other published results using potentials which were optimized to a small family of molecules. The calculated structures for the azabenzenes are very similar to those of Williams and Weller,5 which is not surprising since the two models mainly differ in the use of a full DMA, rather than potential derived charges with explicit lone-pair sites, to represent the electrostatic potential. 4. Discussion The electrostatic forces between polar and aromatic molecules clearly play an important role in determining their crystal structures. The large quantitative, and sometimes qualitative, changes in the calculated structures that occur when the anisotropic multipoles are removed demonstrate clearly that the packing of most of these molecules is sensitive to the detailed model for the electrostatic forces. This implies that different atomic charge models could give very different results. Nevertheless, the demonstration by Williams and Weller5 that the azabenzene structures could not be reproduced by even a potential derived charge model unless additional lone-pair sites were used, shows that no atomic charge model will be capable of calculating reliably the crystal structures of a wide range of molecules. This is because it is incapable of accurately representing the electrostatic potential in the van der Waals contact region. Simpler electrostatic models than the DMA will be adequate for modeling the crystal structures of many molecules, but it is clearly necessary to use an electrostatic model that represents the electrostatic forces accurately for confidence that the scheme will extend to a wide range of molecules and crystal structures. This sensitivity of the calculated crystal structure to the detailed form of the electrostatic model has also been seen for some explosive molecular crystals.17 Although a distributed multipole representation of the electrostatic forces is far less sensitive to the basis set used for the wave function than most atomic charge models, there is still a variation resulting from the change in the quality of the wave function. Our results demonstrate that the strength of the electrostatic forces relative to the other components of the potential, and therefore the quality of the wave function used to derive the electrostatic model, does affect the calculated structure as well as the lattice energy. The scaling factor of 0.9 to allow for the effect of electron correlation on SCF 6-31G** wave functions is a crude approximation, as this ratio varied from 0.78 for NH3 to 0.96 for CH3CN for the molecules considered by Cox and Williams,23 and this factor may be even more approximate for quadrupolar molecules. It would clearly be desirable to use the DMAs of high-quality, correlated wave functions to model the crystal structures of the more sensitive molecules. With the most realistic electrostatic model, most of the lattice parameters for the crystal structures are being reproduced to
J. Phys. Chem., Vol. 100, No. 18, 1996 7359 within a few percent. The calculations use conventional static lattice minimization which in principle calculates the 0 K structure, yet the comparison is made with the room temperature one. This method has previously been justified by the fact that the thermal effects are minor relative to the errors produced by the potential model and that they may be partially absorbed into the fitted repulsion-dispersion parameters. However, crude comparisons with the low-temperature crystal structures which exist for some of the molecules in our selection34 show that the cell edges contract usually on the order of 1-2% on cooling but that this is extremely variable and anisotropic, with individual cell edges contracting by up to 4% (p-benzoquinone) or even expanding by nearly 2% (formamide) on cooling. Hence, even though some thermal effects will have been absorbed in the empirical fitting of the repulsion-dispersion parameters, this cannot adequately account for the temperature effects. Similarly, the approximation of a rigid molecule will be particularly poor for the NH2 groups which are partially sp3 hybridized, and the determination of a mean position for the mobile protons of this functional group by X-ray crystallography is difficult. Thus, a significant proportion of the larger errors noted for many of the molecules with sp3 amino groups can be attributed to the errors in representing this group (Table 4). Thus, it seems likely that a large proportion of the errors remaining in these crystal structure calculations come from the static simulation technique. There is still room for improvement in the model potential, as the failure of our attempts at potential optimization probably indicates that our initial parameters are only about right for such a wide range of molecules within the implicit assumptions of transferability and functional form. The potential could be improved for specific molecules or families, by using the best available models for all components of the potential, where there are sufficient experimental and theoretical data available. 5. Conclusion The potential schemes EST+0.9DMA and FIT+0.9DMA, based on an accurate electrostatic model, provide a transferable approach to deriving a potential that is adequate for reproducing the room temperature crystal structures of many rigid organic molecules within the limits of a static minimization method. The anisotropy of the electrostatic forces sufficiently dominates the molecular packing in many cases, so that a reasonably accurate description provides a potential scheme that is often a major improvement on isotropic atom-atom potentials. Acknowledgment. This work was supported by the EPSRC under Grants GR/J02872 and GR/F76718, with contributions from The Wellcome Foundation Limited. The computer resources were provided by EPSRC under Grant GR/J31865. References and Notes (1) Pertsin, A. J.; Kitaigorodsky, A. I. The Atom-Atom Potential Method; Springer-Verlag: Berlin, 1987. (2) Filippini, G.; Gavezzotti, A. Acta Crystallogr. 1993, B49, 868880. (3) Gavezzotti, A.; Filippini, G., J. Phys. Chem. 1994, 98, 4831-4837. (4) Desiraju, G. R. Crystal Engineering: The design of organic solids; Elsevier: Amsterdam, 1989; p 56. (5) Williams, D. E.; Weller, R. R. J. Am. Chem. Soc. 1983, 105, 41434148. (6) Price, S. L. Mol. Phys. 1987, 62, 45-63. (7) Buckingham, A. D.; Fowler, P. W. Can. J. Chem. 1985, 63, 20182025. (8) Price, S. L.; Stone, A. J. J. Chem. Phys. 1987, 86, 2859-2868. (9) Mitchell, J. B. O.; Nandi, C. L.; Thornton, J. M.; Price, S. L.; Singh, J.; Snarey, M. J. Chem. Soc., Faraday Trans. 1993, 89, 2619-2630.
7360 J. Phys. Chem., Vol. 100, No. 18, 1996 (10) Price, S. L.; Lo Celso, F.; Treichel, J. A.; Goodfellow, J. M.; Umrania, Y. J. Chem. Soc., Faraday Trans. 1993, 89, 3407-3417. (11) Hurst, G. J. B.; Fowler, P. W.; Stone, A. J.; Buckingham, A. D. Int. J. Quantum Chem. 1986, 29, 1223-1239. (12) Berkovitch-Yellin, Z.; Leiserowitz, L. J. Am. Chem. Soc. 1980, 102, 7677-7690. (13) Berkovitch-Yellin, Z.; Leiserowitz, L. J. Am. Chem. Soc. 1982, 104, 4052-4064. (14) Stone, A. J.; Alderton M. Mol. Phys. 1985, 56, 1047-1064. (15) Allen, F. H.; Kennard, O. Chem. Des. Autom. News 1993, 8, 3137. (16) Allen, F. H.; Kennard, O.; Watson, D. G.; Brammer, L.; Orpen, A. G.; Taylor, R. J. Chem. Soc., Perkin Trans. 2 1987, S1-S19. (17) Ritchie, J. P.; Kober, E. M.; Copenhaver, A. S. Proceedings of the 10th International Symposium on Detonation, Office of Naval Research, 1993. (18) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213222. (19) CADPAC5: The Cambridge Analytic DeriVatiVes Package Issue 5, A suite of quantum chemistry programs developed by RD Amos with contributions from I. L. Alberts, J. S. Andrews, S. M. Colwell, N. C. Handy, D. Jayatilaka, P. J. Knowles, R. Kobayashi, N. Koga, K. E. Laidig, P. E. Maslen, C. W. Murray, J. E. Rice, J. Sanz, E. D. Simandiras, A. J. Stone, and M.-D. Su, Cambridge, U.K., 1992. (20) Price, S. L.; Stone, A. J.; Alderton, M. Mol. Phys. 1984, 52, 9871001. (21) Smith, W. The program MDMULP, CCP5 program library; Daresbury Laboratory: Warrington, England, 1982. (22) Ryan, M. D. Modeling the Hydrogen Bond. ACS Symp. Ser. 1994, 569, 36-59. (23) Cox, S. R.; Williams, D. E. J. Comput. Chem. 1981, 2, 304-323. (24) Price, S. L.; Andrews, J. S.; Murray, C. W.; Amos, R. D. J. Am. Chem. Soc. 1992, 114, 8268-8276.
Coombes et al. (25) Price, S. L. Philos. Mag. 1996, 73, 95-106. (26) Williams, D. E.; Cox, S. R. Acta Crystallogr. 1984, B40, 404417. (27) Cox, S. R.; Hsu, L.-Y.; Williams, D. E. Acta Crystallogr. 1981, A37, 293-301. (28) Mitchell, J. B. O.; Price, S. L. J. Comput. Chem. 1990, 11, 12171233. (29) Stone, A. J.; Tong, C.-S. J. Comput. Chem. 1994, 15, 1377-1392. (30) Williams, D. E.; Starr, T. L. Comput. Chem. 1977, 1, 173-177. (31) Willock, D. J.; Price, S. L.; Leslie, M.; Catlow, C. R. A. J. Comput. Chem. 1995, 16, 628-647. (32) To illustrate this problem, a revision of the sublimation of cytosine by 21 kJ mol-1 came to our attention during the course of our work, bringing the experiment into agreement with our predictions. We have also not included a lattice energy of -57.3 kJ mol-1 for TRPHAM in Table 3, as this seems unlikely to refer to complete separation of the molecules. (33) Price, S. L.; Stone, A. J. Mol. Phys. 1984, 51, 569-583. (34) Price, S. L. In Computer Modelling in Inorganic Crystallography; Catlow, C. R. A., Ed.; Academic Press: London, in press. (35) Hagler, A. T.; Lifson, S.; Dauber, P. J. Am. Chem. Soc. 1979, 101, 5122-5130. (36) Filippini, G.; Gavezzotti, A. Chem. Phys. Lett. 1994, 231, 86-92. (37) Chickos, J. S. In Molecular Structure and Energetics; Liebman, J. F., Greenberg, A., Eds.;VCH: New York, 1987; Vol. 2, p 67. See also personal communication of updated tables. (38) Nabavian, M.; Sabbah, R.; Chestel. R.; Laffite, M. J. Chim. Phys. Phys.-Chim. Biol. 1977, 74, 115-126. (39) Bondi, A. J. Chem. Eng. Data 1963, 8, 371-381.
JP960333B
