A novel definition of a molecule in a crystal

~r-m--'g"V,~m

21 March 1997

CHEMICAL PHYSICS LETTERS

ELSEVIER

Chemical Physics Letters 267 (1997) 215-220

A novel definition of a molecule in a crystal Mark A. Spackman *, Patrick G. Byrom Department of Chemistry, University of New England, Armidale NSW 2351, Australia

Received 17 December 1996; in final form 17 January 1997

Abstract A new method for dividing a crystalline electron distribution into molecular fragments is proposed, based on Hirshfeld's partitioning scheme. Unlike other approaches, the method partitions the crystal into smooth molecular volumes as well as intermolecular voids of tow electron density. To compare the new method with several other schemes which subdivide a crystal into molecules, numerical integration is performed on two model electron densities (one representing a superposition of isolated molecules, the other interacting molecules) for ice VIII, formamide and urea. The new scheme is simple to apply, aesthetically appealing, and offers some promise in routine partitioning of crystalline electron densities or in computer graphics to provide additional insight into molecular packing in crystals.

1. Introduction The definition of a molecule in a condensed phase is an elusive goal, yet one which continues to attract attention because the recognition of distinct entities in molecular liquids and crystals is a fundamental concept in chemistry. For example, the molecular dipole moment of water is an important parameter in the realistic simulation of its dielectric properties [ 1], and its estimation in the condensed phase using various models has been of interest for several decades [2-6]. Our interest lies in the extraction of molecular properties from accurate x-ray diffraction data [7,8] for crystals, and for that purpose the definition of a molecule as a fragment of the total crystalline electron distribution is unavoidable, yet acknowledged to be arbitrary [9-12]. We have been exploring various schemes which partition the crystalline electron density into molecular fragments, and * Corresponding author.

here we report a scheme based on Hirshfeld's stockholder partitioning [13]. The scheme is simple to apply, aesthetically appealing, and offers some promise as a pedagogical tool in discussions of molecular packing in crystals. Our discussion and presentation will focus on electron distributions derived from modelled x-ray diffraction data [10,11,14], although the conclusions obtained are of a more general nature. To place the new partitioning scheme in sufficient perspective we compare and contrast it with other schemes which have been used for this purpose. Those of direct relevance to our main goal, including the new scheme, are outlined: (i) Generalized Wigner-Seitz (WS) partitioning is based on the original concept of an atomic cell (as applied to sodium [15]), generalized to apply to molecules and used by Coppens and co-workers to estimate charges on molecular ions [16,17] and dipole and quadrupole moments for a number of molecules [18-20], all directly from x-ray diffraction data. The

0009-2614/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved. PH S0009-261 4(97)00100-0

M.A. Spackman, P.G. Byrom / Chemical Physics Letters 267 (1997) 215-220

216

molecular surface is defined by a set of conditions of the form r A • UA._._.~B_ RA

] --

rB" UAB

I

RB

1'

where r A and r B are vectors from atoms A and B in two different molecules respectively, UAB is the unit vector from nucleus A to nucleus B, and R A and R B are appropriate van der Waals radii of the atoms. These conditions yield a molecular surface as the union of planar segments. A simpler form has often been implemented, eliminating the projection on the unit vector:

and these conditions yield a molecular surface as a union of hyperbolic sections. The molecular WS surface is thus readily defined, but its precise nature depends strongly on the atomic radii chosen, and the surfaces generated are at best piecewise smooth; an example for formamide has been reported [19], and one for urea (using the simpler criteria, see details below) is presented in Fig. lb. The precise surface also depends on whether the criterion is applied pairwise to atoms in the molecule and the closest atoms in neighbouring molecules, or whether the test is made over aggregates of atoms. (ii) The quantum theory of atoms in molecules (QTAM) has been described in detail by Bader [21],

a

b

and is being employed increasingly in research in chemical structure and bonding. It has only recently been applied to the extraction of molecules in crystals, in particular urea [22], ice [2] and hydrogen cyanide [23]. For this purpose the molecule is the union of its atomic basins, each basin enclosed by the surface defined by

Vp(r) n(r) = o, where n(r) is a unit vector normal to the surface at r and Vp(r) the gradient of the electron density in the crystal, p(r). The definition of atomic basins has a rigorous quantum mechanical foundation, and their union to form molecules yields molecular surfaces in crystals which are piecewise smooth; an example for urea is given in Fig. lc. Unlike the other schemes considered here, Bader's scheme depends intimately on the crystalline electron density p(r); comparison of molecular properties obtained from different p(r) will involve different molecular surfaces. As discussed in recent work on urea [22], the numerical integration of molecular properties using this partitioning scheme is difficult and time-consuming because of the practical problem of accurately defining the molecular surface in crystals. Preliminary application to experimental electron densities has recently been reported [24] but the algorithm used is presently very time-consuming. (iii) Hirshfeld partitioning ( H ) is an extension of Hirshfeld's stockholder concept [13] which divides

c

Fig. 1. Sections in the molecular plane of different molecular surfaces for urea: (a) the shaded area is inside the HS surface defined by w(r) = 0.5, the two other contours of w(r) shown are 0.1 and 0.9; (b) the shaded area is inside the WS surface; (c) an example of the Bader surface from an ab initio crystal Hartree-Fock calculation on urea [22] (the atomic and molecular surfaces are defined by the locus of gradient paths emanating from each of the atomic nuclei). All maps are 8 ,~ square.

M.A. Spackman, P.G. Byrom/ Chemical Physics Letters 267 (1997) 215-220

the electron density of a molecule into continuous atomic fragments. The concept was generalized to extract overlapping, continuous molecular fragments from experimental electron densities in crystals by defining a molecular weight function: w(r) =

E A ~ molecule

PA(r)/

~.,

PA(r)

217

1 0.5 0 -0.5 -1

A ~ crystal -4

= Ppromolecule(r)/Pprocrystal(r).

Here pA(r) is a spherically-averaged atomic electron density function [25] centred on nucleus A, and the promolecule and procrystal are appropriate sums over the atoms belonging to a single molecule and the crystal, respectively, w(r) satisfies the condition 0 < w(r) < l, and molecular properties are obtained by integration over the weighted electron density, w(r)p(r). Using this form of partitioning, molecular properties have been reported for several compounds, either directly from the experimental structure factors [18-20] or by partitioning of the multipole-refined electron density in direct space [26-29]. The most significant difference between this scheme and the WS and QTAM schemes is the absence of a molecular surface, which simplifies the integration for determining molecular properties. (iv) The new scheme, based on Hirshfeld's partitioning, is denoted Hirshfeld partitioning with a surface (HS) as it uses the same weight function w(r) in a slightly different manner. Instead of weighting the crystalline electron density by w(r) we define a molecule as that region where w(r)>10.5 (i.e. the promolecule contribution exceeds that from all neighbouring molecules). Within this region all of the electron density is counted for purposes of integration; outside the region none is counted. Fig. l a presents three contours of w(r) in the molecular plane for urea, and the region inside the 0.5 contour - the molecular region - is shaded. The relationship to the normal Hirshfeld partitioning is readily seen, and the broad similarity to the piecewise smooth WS surface (Fig. lb) is obvious. Fig. 2 depicts w(r) for urea as a relief map, in the same plane as Fig. l, and here we see rather dramatically the nature of the molecular weight function w(r): it is extremely flat over most of the molecule, dropping rapidly beyond approximately 0.5 ,~ from any of the nuclei. The cutoff value of w(r)= 0.5 used to define the HS

-3 -2 -1 0

-4

-3

-2

-1

0

1

2

3

4

Fig. 2. Relief and contour maps of the Hirshfeld weight function w(r) in the molecular plane of urea (same data as Figure la), which demonstrate the flat nature of w(r) in the vicinity of the atoms in the molecule (where w(r) is very close to 1.0), and the reasonably sharp decrease with distance from the molecule.

surface is not arbitrary; it guarantees maximum proximity of neighbouring molecular volumes, yet prevents overlap. A cutoff of less than 0.5 would yield more space between adjacent molecules, while a cutoff greater than 0.5 would allow molecular volumes to overlap. The HS molecular surface is everywhere smooth and its definition relatively trivial, especially for numerical integration and computer graphics. The most conspicuous difference between this scheme and the other three above (and in part its novelty) is that it partitions the crystal into molecular regions as well as intermolecular 'voids'. We show below that the interstices, where no single molecule dominates Pprocrystal(r),are almost devoid of electron density.

2. Results of numerical integration

As the QTAM scheme has not yet been applied routinely to the type of electron densities we are studying, we have used the WS, H and HS schemes to extract molecular fragments and their dipole moments from crystalline electron densities for urea, formamide and ice VIII. For these purposes we constructed two different model electron densities for each system (molecules and crystal) by refine-

218

M.A. Spackman, P.G. Byrom / Chemical Physics Letters 267 (1997) 215-220

ment of a pseudoatom multipole model via least squares fitting to structure factors generated from ab initio electron densities for a superposition of isolated molecules in the crystal (molecules) and a crystal Hartree-Fock calculation (crystal). These crystalline electron distributions consist of a superposition of realistic analytical (smooth and continuous) molecular electron density functions with known dipole moments, and the differences between the molecules and crystal dipole moments are significant and reflect the polarization of the electron density in these hydrogen-bonded systems. Our interest lies in the extent to which the various space-partitioning schemes retrieve the molecular dipole moments and the difference between the molecules and crystal values. Details of the model electron densities are unimportant for our present purpose, and will be published elsewhere [30]. Radii of 1.14 ,~ (H), 1.68 .~ (C), 1.54 ,~ (N), and 1.46 ,~ (O) [31] were used to define the WS surface, and the H and HS schemes employed sphericallyaveraged Hartree-Fock limit electron densities [25]. Numerical integration utilized a three-dimensional cartesian grid with resolutions of (in .~) 0.050 × 0.050 × 0.050 (urea and formamide) and 0.033 × 0.033 × 0.029 (ice VIII). Although no formal error

analysis was performed in these preliminary studies, the number of electrons in the molecule, Nel (Table 1), is one measure of the accuracy of the integration, and in all instances the integral is within 0.5% of the expected value. We estimate integration errors of approximately 0.05 D in each of the results in Table 1, and suspect this is rather conservative. Main outcomes of the integration can be summarized pointwise: (a) all three schemes generally underestimate the true values; (b) WS and HS results are similar, but the HS result is most often closest to the true value (typically 13% low) while the H result is always furthest away (typically 33% low); (c) in spite of (a) all three schemes yield similar estimates of the dipole moment enhancement (WS and HS results are almost identical), although none of these is very accurate; (d) within the accuracy of our integration the intermolecular voids in the HS scheme contain no significant electron density. Although only preliminary, these results are of considerable significance to efforts to extract molecular properties from x-ray diffraction data, and shed light on the varying results obtained from using

Table 1 Molecular dipole moment magnitudes (Debye) from numerical integration over molecular volumes defmed by generalised Wigner-Seitz (WS), Hirshfeld (H) and Hirshfeld surface (HS) partitioning schemes. Nej is the number of electrons in the molecular volume, and dipole moments are reported from electron densities representing a superposition of isolated molecules (molecules) and a crystal Hartree-Fock result (crystal) Nej

Molecules

Crystal

Crystal - molecules

WS H HS true value

32.12 32.13 32.08 32.00

3.13 2.26 3.28 4.38

4.31 3.37 4.47 5.16

1.18 1.11 1.19 0.78

formamide WS H HS true value

23.92 23.89 23.92 24.00

3.50 2.77 3.11 3.66

4.52 3.66 4.13 4.41

1.02 0.89 1.02 0.75

ice VIII WS H HS true value

10.00 9.98 9.95 10.00

1.43 1.05 1.53 1.64

1.78 1.36 1.88 2.14

0.35 0.31 0.35 0.50

urea

M.A. Spackman, P.G. Byrom / Chemical Physics Letters 267 (1997) 215-220

219

Fig. 3. Examples of molecules in crystals defined by the HS scheme, left to right: urea, formamide and ice VIII. Top: a single molecule (shaded) superimposed on a contour map of the molecules electron density used for integration purposes (electron density contours are at values of 2 - " e A - 3 for n = - 5 to n = 4). Bottom: packing diagrams showing the close approach of nearest neighbour molecules and the intermolecular regions belonging to no molecules.

different space-partitioning schemes reported in the literature [7]. In particular, systematically low estimates provided by the Hirshfeld scheme (H) have been noted, and these arise from the tendency of the weighting function w(r) to bias the molecular property towards the promolecule result (e.g. zero dipole moment). The normal Hirshfeld scheme would seem to be a poor choice for the extraction of molecules from crystals.

3. Packing of molecular volumes in the crystal The potential for use of the present HS partitioning scheme in molecular packing discussions is highlighted by Fig. 3. For each system in the figure the HS surface for an individual molecule passes along or close to the regions of lowest electron density between neighbouring molecules, in this manner closely approximating Bader's scheme, for which the

Fig. 4. Three-dimensional isosurfaces (w(r) = 0.5) of HS-defined molecules in the crystals of (left to right) urea, formamide and ice VIII. These images were created using the isosurface facility of the NCSA XDataSlice program [32].

220

M.A. Spackman, P.G. Byrom / Chemical Physics Letters 267 (1997) 215-220

molecular surface in any plane necessarily coincides with the valleys of the electron distribution. The molecular packing diagrams confirm the close approach of the HS molecular fragments, especially along hydrogen bond directions. These diagrams offer additional insight into the packing of molecules in crystals, and afford a hitherto unseen picture of molecular shape in a crystalline environment. Fig. 4 provides a shaded rendering of individual molecules extracted from their crystal lattice. The HS-defined molecules are rather beautiful entities, remarkably smooth but faceted, reflecting the proximity of near neighbours. It should prove extremely interesting to encode a property such as the electrostatic potential on these surfaces, or to use the molecular volumes in routine partitioning of crystalline electron densities. Although in this preliminary communication we have only explored the application of the novel Hirshfeld surface scheme to three molecular crystals, we believe this simple definition of a molecule in a crystal will lend itself to a wide variety of future applications. With reference to molecular properties, we plan to explore its similarities with, and differences from, QTAM partitioning, since the most significant differences between the two schemes appear to be in regions of very low electron density. Finally, we note that the weighting function does not necessarily need to b e constructed from spherical atoms, but rather any spherical shape function which sensibly reflects relative atomic sizes should suffice, and give a subtly different molecular surface. Acknowledgements This work has been supported by the Australian Research Council. We are indebted to Dr C. Gatti (CNR, Milano) for providing us with Fig. 1c. References [1] M. Sprik, J. Chem. Phys. 95 (1991) 6762. [2] C. Gani, B. Silvi and F. Coionna, Chem. Phys. Lett. 247 (1995). [3] P. Barnes, J.L. Finney, J.D. Nicholas and J.E. Quinn, Nature 282 (1979) 459,

[4] D.N. Bernardo, Y. Ding, K. Krogh-Jespersen and R.M. Levy, J. Phys Chem. 98 (1994) 4180. [5] C.A. Coulson and D. Eisenberg, Proc. R. Soc. London, A 291 (1966) 445. [6] K. Laasonen, M. Sprik, M. Parrinello and R. Car, J. Chem. Phys. 99 (1993) 9080. [7] M.A. Spackman, Chem. Rev. 92 (1992) 1769. [8] M.A. Spackman and P.G. Byrom, Acta Crystallogr., Sect. B 52 (1996) 1023. [9] R.F. Stewart, Chem. Phys. Lett. 65 (1979) 335. [10] F.L. Hirshfeld, Cryst. Rev. 2 (1991) 169. [ 11] F.L. Hirshfeld, in: Accurate Molecular Structures. Their Determination and Importance, eds. A. Domenicano and I. Hargittai (IUCr, Oxford Univ. Press, Oxford, 1992) p. 237. [12] F.L. Hirshfeld, J. Mol. Struct. 130 (1985) 125. [13] F.L. Hirshfeld, Theor. Chim. Acta 44 (1977) 129. [14] G.A. Jeffrey and J.F. Piniella, The application of charge density research to chemistry and drug design (Plenum Press, New York, 1991). [15] E. Wigner and F. Seitz, Phys. Rev. 43 (1933) 804. [16] P. Coppens, Phys. Rev. Lett. 35 (1975) 98. [17] P. Coppens and T.N. Guru-Row, Ann. N.Y. Acad. Sci. 313 (1978) 244. [18] P. Coppens, G. Moss and N.K. Hansen, in: Computing in Crystallography, eds. R. Diamond, S. Ramaseshan and K. Venkatesan (Indian Academy of Sciences, Bangalore, 1980) p. 16.01. [19] G. Moss and P. Coppens, Chem. Phys. Lett. 75 (1980) 298. [20] G. Moss, in: Electron Distributions and the Chemical Bond, eds. P. Coppens and M.B. Hall (Plenum Press, New York, 1982) p. 383. [21] R.F.W. Bader, Atoms in Molecules - a Quantum Theory (Oxford Univ. Press, Oxford, 1990). [22] C. Gatti, V.R. Saunders and C. Roetti, J. Chem. Phys. 101 (1994) 10686. [23] J.A. Platts and S.T. Howard, J. Chem. Phys. 105 (1996) 4668. [24] C. Flenshurg, D. Madsen, S. Larsen and R.F. Stewart, Acta Crystallogr., Sect. A 52, Supplement (1996) C. [25] E. Clementi and C. Roetti, At. Data Nucl. Data Tables 14 (1974) 177. [26] F.L. Hirshfeld and H. Hope, Acta Crystallogr., Sect. B 36 (1980) 406. [27] F.L. Hirshfeld, Acta Crystallogr., Sect. B 40 (1984) 484. [28] F. Baert, P. Coppens, E.D. Stevens and L. Devos, Acta Crystallogr., Sect. A 38 (1982) 143. [29] M. Eisenstein, Aeta Crystallogr., Sect. B 44 (1988) 412. [30] M. Alfredsson, P.G. Byrom, K.G. Hermansson and M.A. Spackman, in preparation (1996). [31] M.A. Spackman, J. Chem. Phys. 85 (1986) 6579. [32] XDataSlice 2.2 (Software Development Group, National Centre for Supercomputing Applications, University of Illinois at Urbana-Champaign, 1993).