Critical Evaluation of C−H···X Hydrogen Bonding in the Crystalline

Kikuko Hayamizu , Seiji Tsuzuki and Shiro Seki. The Journal of Physical ... Yumi N. Imai, Yoshihisa Inoue, and Yoshio Yamamoto. Journal of Medicinal ...
0 downloads 0 Views 1MB Size
CRYSTAL GROWTH & DESIGN 2003 VOL. 3, NO. 5 643-661

Articles Critical Evaluation of C-H‚‚‚X Hydrogen Bonding in the Crystalline State Jan-Albert van den Berg* and Kenneth R. Seddon* The QUILL Centre, The Queen’s University of Belfast, United Kingdom, BT9 5AG Received May 21, 2003;

Revised Manuscript Received June 30, 2003

ABSTRACT: Contacts of the type C-H‚‚‚X where X ) F, Cl, Br, I, O, S, or N were analyzed statistically using structures available in the Cambridge Structural Database. For this analysis, the cone correction method was extended to an isotropic density correction. This method corrects for both angular and distance effects by normalizing isotropic contact densities from the experimental distributions. It was possible to show the widespread occurrence of C-H‚‚‚X interactions in the crystalline state using this method. Distances of highest incidence, at which the largest number of contacts were sampled after normalization of isotropic distribution, were determined for each interaction type. These distances were compared to the classically accepted van der Waals cutoff radii, often used to distinguish between hydrogen-bonded interactions and van der Waals contacts. Generally, distances of highest incidence did not occur at the sum of van der Waals cutoff radii but occurred either below or above that distance depending on the specific interaction analyzed. This methodology, while of course not yielding quantitative data on hydrogen bond strength, allowed for a qualitative ordering of the interactions studied, in terms of expected interaction strengths. Introduction The phenomenon of “the hydrogen bonding interaction” was described by Pauling1 as “...an atom of hydrogen [...] attracted by rather strong forces to two atoms, instead of one, so that it may be considered [...] as a bond between them.” This situation, because it is mainly attributed to electrostatic forces, contrasts with the interactions between species due to dispersion forces, normally referred to as a van der Waals interaction. However, the limits of hydrogen bonding in general have been an area of continuous debate, due to difficulties in experimentally verifying all possible hydrogenbonded interactions between individual interacting species. The chemical bond is defined2 by its ability to stabilize an aggregate of atoms, and so becomes an independent molecular entity, and intermolecular bonds have required a very careful approach to bond definition due to their varied strengths and interacting distances. 3 Lately, the number of articles that describe hydrogen bonds which are much weaker than the classically accepted hydrogen bonds has increased considerably, and the distinction between the “van der Waals contact” and the “hydrogen bond” has been a topic of serious reevaluation. Steiner (who has recently written a com* To whom correspondence should be addressed. j.vandenberg@ qub.ac.uk; [email protected].

prehensive review4 on the latest developments of hydrogen bonding in the crystalline state) and Desiraju5 proposed that the distinction between weak hydrogen bonds and van der Waals contacts, specifically as observed from distances in the crystalline state, varies between those expected for each extremesit is misleading to accept that all hydrogen bonds (or van der Waals contacts for that matter) are alike. However, this approach is not generally accepted, and indeed Cotton et al.6 stated: It is clear that the field is getting muddier and muddier as the definition of a hydrogen bond is relaxed...we strongly disagree with the newer and more relaxed definitions that do not distinguish between a ‘hydrogen bond’ and what is nothing more than a classical van der Waals interaction. When two atoms, ions, molecules, or aggregates approach each other, the equilibrium distance between them (which is the sum of their individual van der Waals radii, rXw + rYw, where X and Y are the atoms in closest contact) is a result of the balance between the attractive and repulsive forces between the units. These forces are normally described by the Lennard-Jones potential,7 attractive at long distances as r-6 and repulsive at short distances as r-12 (see Figure 1). Bondi8

10.1021/cg034083h CCC: $25.00 © 2003 American Chemical Society Published on Web 07/19/2003

644

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

van den Berg and Seddon Table 1. Summary of Bondi’s van der Waals Radiia atom, X

aggregate

rXw/Å

atom, X

aggregate

rXw/Å

F Cl Br I H O O

C-F C-Cl C-Br C-I C-H C-O-C CdO

1.50 1.76 1.85 1.96 1.20 1.52 1.50

S N N N P C C

C-S-C NH3 CdN-C CtN P4 or PH3 Caliph Carom

1.83 1.55 1.55 1.60 1.80 1.70 1.77

a

Ref 8.

Scheme 1. Summation of the van der Waals Radii of Interacting Atoms Used to Determine Whether a Hydrogen Bond Occurred Figure 1. The Lennard-Jones potential curve7 describes the balance between attractive and repulsive forces between molecular aggregates.

Figure 2. The van der Waals radii (rXw and rYw) for atoms in a heteroatomic diatomic molecule. These radii are used in volume calculations to ensure transferability.

developed a scheme, now used universally, to determine undefined van der Waals radii from empirical crystallographic data, with the implicit (but demonstrably incorrect) assumption that the values of rw are both constant and transferable for any aggregate situation. By using radii determined from crystallographic data to predict van der Waals volumes (as shown in Figure 2), it is however possible to determine upper and lower limits between which the van der Waals radii must lie. In addition to the fallacy that the radii so obtained be transferable between aggregates, Bondi made another dubious assumption, albeit minor: a spherical model was assumed for the atoms without any correction to 0 K (a rather laborious process for larger aggregates and not necessary if an upper limit on rw is to be found). These values from Bondi have been invariably applied to crystallographic data as a means of “determining” the existence or nonexistence of a possible hydrogen bond. Because the hydrogen bond is thermodynamically stronger than the van der Waals interaction, it appears logical to assume that hydrogen bonds will be shorter in distance, an incorrect assumption. It is thus argued that an easy to apply cutoff criterion for establishing the existence of a hydrogen bond is simply to add the van der Waals radii of the interacting atoms, and check whether the experimentally observed interaction occurs over a shorter distance than that of the calculated one (see Scheme 1). If it is the case, then it follows that the interaction is hydrogen bonded. The most common radii used to establish this criterion are listed in Table 1. The application of hydrogen bonds in protein chemistry and DNA structure elucidation is well-known, and its importance in the field of biochemistry cannot be

underestimated.9-11 However, as a consequence of the distorted message from van der Waals radius considerations, an implicit assumption in the application of hydrogen bonds to theoretical calculations denied the existence of the C-H‚‚‚X hydrogen bonds without differentiation, implying that these interactions do not have any effect on the structural properties of the aggregates in study.12,13 Early research14 showed that C-H groups might be involved in hydrogen bonding interactions, and later the C-H‚‚‚O interaction15,16 was clearly illustrated17 to impart structural changes to aggregates in which it occurs. This immediately pointed to the need for a change in the philosophy adopted for these interactions: why should only some C-H‚‚‚X interactions be included as hydrogen bonds while others are ignored? Is it possible that environmental changes due to different atoms involved (e.g., C-H‚‚‚Cl- vs. C-H‚‚‚Cl-M vs. C-H‚‚‚Cl-R) do not alter the properties of the interaction concerned? Here it is proposed that C-H‚‚‚X hydrogen bonds are real, that they do play an important part in the crystallization of myriads of organic aggregates, and will be an important factor in the prediction of new and novel crystalline aggregates (predictable supramolecular synthons).18-20 This is not to say that all C-H‚‚‚X contacts are hydrogen bonds, only that many interactions are ignored when the classical model is used. Hydrogen bonds have been characterized by their directional15 and electrostatic21 (as r-2, and therefore occurring at shorter distances) properties. In this paper, focus will be on the directional property of hydrogen bonds: the interaction between the available electron density on the donor atom (classically described as a lone pair for such atoms) and vacant orbitals on the acceptor atoms causes a tendency to align D-H‚‚‚A (with D the donor and A the acceptor) in a near linear fashion. Due to other intermolecular forces at play during the process of crystallization, it is rarely the case that an exactly linear interaction is observed, thus the need to account for angular distribution effects (normally presented as histograms of distance vs. angle) and the development of the cone correction method. The Cambridge Structural Database22,23 (CSD) has enabled various researchers12,14,16,24-27 to study the effect of hydrogen bonds, in crystalline aggregates, across a wide range of structures. Kroon and Kanters28 developed the well-known cone correction to histograms obtained, to normalize angular distribution effects.

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Figure 3. Analysis of the parameters requested during the contact search; (a) schematic of the actual parameters (r,φ); (b) the cone that emanates from the hydrogen atom leading to the cone correction method; (c) the area element which, after rotation about the D-H axis yields a ring-like volume element, Vij, the actual volume sampled during contact counting. The hydrogen atom is always placed at the origin.

Using the cone correction method on angular histograms clearly showed that, while typical van der Waals interactions, e.g., R3C-H‚‚‚H-CR3, do not have any angular preference; stronger hydrogen bonded interactions, e.g., RO-H‚‚‚Cl-, mainly aligned in a linear orientation. It is the purpose of this paper to study the effectiveness of the van der Waals criterion as an approach to prove or disprove the existence of C-H‚‚‚X hydrogen bonds in the crystalline phase, and so analyze the validity of Steiner and Desiraju’s proposition5 by looking at various hydrogen-bonded contacts from a statistical point of view. The cone correction method is extended to the isotropic density correction (IDC), which factors out isotropic nonbonded atomic distribution at distances far away from the hydrogen atom (thus, normalizing both angular and distance distribution). It is demonstrated that using only empirical distance and angle relationships between a large number of interacting aggregates, a clear distinction between van der Waals interactions (which are by definition not directional) and hydrogen bonded interactions exists. This paper builds on a preliminary communication,12 and now includes an improved analytic approach and a much wider empirical data set. Experimental Section Contact Searching. The CSD22,23 version 5.22 (245392 entries) was used to search for various interaction classes. All searches required that the R-factor be less than 0.05, and that there be no disorder or errors in the crystal structures, as defined by QUEST.29 For simplicity’s sake, these classes were divided into various interaction types. An example of an interaction class would be C-H‚‚‚X-, if X- equals a halide anion, whereas a specific interaction type would be C-H‚‚‚Cl-. The program QUEST,29 part of the CSD, was used to construct each interaction search since it offers a user-friendly interface. The contact parameters (r,φ) searched for in each interaction search and ultimately stored in a table file (.tab) are shown in Figure 3a. Notice that the distance between the donor atom and the hydrogen atom is not stored since hydrogen atom positions cannot be accurately determined using X-ray techniques alone (often the C-H distance is arbitrarily fixed). However, this does not affect our study because it is the directionality of these interactions that is important in distinguishing whether an interaction was hydrogen bonded or not. The program VISTA22 can usually be used to analyze the list of data pairs obtained. Since the version of VISTA that was available only allowed for a maximum of 105 contacts to be analyzed (and the number of interactions obtained for a specific type can be much larger), a small program, ISOCOR,

Crystal Growth & Design, Vol. 3, No. 5, 2003 645

Figure 4. Included and excluded hydrogen bonded contacts. The example of a 1,4-butanediol dimer is used merely to serve as an illustration (it was not taken from the CSD). All the red contacts would be excluded (since they are intramolecular and separated by less than six bonds) while the blue contacts would be included. was written in C++ to count the contacts in the .tab file. The program (which is available from the authors) outputs an uncorrected empirical distribution and corrected distributions (refer to the methods of analysis section for a detailed description). The program, KyPlot,30 was used to visualize the matrices (as obtained from ISOCOR) as rainbow contour maps; Excel can be used, but it is much more restrictive in its visualization options. Atoms that were in close proximity due to geometric constraint were excluded from the search (as shown in Figure 4). The search included all intermolecular interactions and excluded all intramolecular interactions separated by less than six bonds. Furthermore, the search was restricted to r/Å ∈ [1,7]; φ/° ∈ [90,180]; and φ/° ∈ [0,360] (the symbols used are shown in Figure 3b for clarity), to restrict the output file size and because contacts outside the r and φ ranges are unlikely to be hydrogen bonds. Methods of Analysis. (a) Data Correction. As described in the Experimental Section, the CSD is searched for all contacts of a certain interaction type (e.g., C-H‚‚‚Cl-C), producing a list of all distance and angle (r,φ) pairs. These pairs were counted by separating them into a grid of distances ∆r and angles ∆φ, so generating an empirical distribution matrix, N, of distances (rows) vs. angles (columns). An element of this matrix will be symbolized by nij, i being the distance index and j the angle index. The cone correction method factors out the effect of the widening radius at the base of the cone shown in Figure 3b. In mathematical terms, if lcirc is the circumference of the circle describing the base of the cone shown in Figure 3b (all equations are solved in radians, but for simplicity sake all the illustrations will be shown in degrees):

lcirc ) 2πr sin(π - φ) ) 2πr sin(φ) ∝ sin(φ)

(1)

Since the probability of finding the acceptor atom on the perimeter of the cone is proportional to lcirc (keeping r constant), we can factor it out of the empirical distribution giving a corrected but not yet normalized distribution. The “cc” above the n indicates “cone-corrected but not normalized”:

ncc ij )

nij sin(φj)

(2)

Notice that the change in the distribution is only dependent on the angle, as indexed by j, occurring in the sine term. Since the distribution thus obtained is a matrix, it must be plotted on three dimensions: angle, distance, and population (such a histogram would be a series of bars looking like a cityscape from above). After the correction has been applied, the adjusted volume of the corrected distribution is not equal to the original, and thus needs to be normalized. A weight factor, W, is calculated to ensure that the volumes of the uncorrected, Vinit,

646 Crystal Growth & Design, Vol. 3, No. 5, 2003

van den Berg and Seddon

and final distributions, Vfin, are equal:

density, Fjiso, can be calculated directly:

∑n , if V

Fjiso )

Since Vinit )

ij

corr

)

i,j

let W )

Vinit Vcorr

∑n

cc ij ,

i,j

cc then nfin ij ) W nij w Vfin ) WVcorr ) Vinit

In the above expression, nfin is the normalized, coneij corrected distribution. It is important to realize that this method can only be applied to a large number of interactions, and has no meaning forsindeed, it cannot even be applied tosspecific contacts. The shape of the actual volume element (after rotation about the D-H axis), dV ) dr dφ dφ, within which interacting contacts for each element of the uncorrected distribution matrix were counted, is illustrated in Figure 3c. To determine the noninteracting contact density within that volume element, the Jacobian31 for the integral change of variable from spherical space coordinates to Euclidean coordinates can be used. Since coordinate conversion occurs via the relationships in eq 3, the Jacobian determinant is eq 5.

x ) r sin(φ) cos(φ) y ) r sin(φ) sin(φ)

(3)

| |

z ) r cos(φ)

∂x ∂r ∂y |Jxyz (r,φ,φ)| ) ∂r ∂z ∂r

|

∂x ∂φ ∂y ∂φ ∂z ∂φ

∂x ∂φ ∂y ∂φ ∂z ∂φ

(4)

sin(φ) cos(φ) r cos (φ) cos(φ) -r sin(φ) sin(φ) r sin(φ) cos(φ) cos(φ) -r sin(φ) 0

) sin(φ) sin(φ) r cos(φ) sin(φ) ) r2 sin(φ)

|

(5)

From eq 5 it therefore follows that

dV ) |Jxyz(r,φ,φ)|dr dφ dφ 2

) r sin(φ)dr dφ dφ

(6)

Since the volume is obtained by rotation (dφ) about the D-H axis (a full rotation is 0 f 2π radians):

dV )



2π 2

0

r sin(φ) dr dφ dφ

) r2 sin(φ) dr dφ





0



(7)

) 2πr2 sin(φ) dr dφ The observant reader might have noticed the recurrence of the sine term throughout the calculations since the change in probability distribution is equivalent to the change in density distribution; the cone correction method is a simplification of the isotropic density correction! However, the r2 term that occurs describes the accumulation of noninteracting isotropic contacts far away from the hydrogen atom. Since the total number of contacts found, Ntot, is known (it is the sum of contacts in the uncorrected distribution matrix), and it is possible to calculate the total volume within which contacts were sampled (a hollow hemisphere surrounding the hydrogen atom since φ/° ∈ [90,180]), Vtot, the average isotropic

Ntot Ntot ) Vtot 2 3 - r3min) π(r 3 max

(8)

The distance range (rmin to rmax), within which contacts were sampled, is specified during the counting process. Thus, the number of expected isotropic contacts, niso ij , within a specific volume element, Vij, can be calculated analytically:

niso j isoVij ) 2πFjiso ij ) F

∫ ∫ ri+1

ri

φj+1

φj

r2 sin(φ) dr dφ

2 3 ) πFjiso(ri+1 - r3i )[cos(φj) - cos(φj+1)] 3

(9)

The integrals in eq 9 describe the small volume element shown in Figure 3c after rotation about the D-H axis, and thus:

niso j isoVij ) ij ) F

Ntot Vij V ) Ntot Vtot ij Vtot

(10)

Thus, finding the relation between the ring-like volume element and the total volume, and multiplying the result with the total number of contacts sampled, determines the isotropic density. Figure 5 illustrates rainbow contour maps (RCMs) derived for a perfect noninteracting contact before and after application of the cone correction method. To derive Figure 5a, eq 9 was calculated over a distance and angle range shown in Figure 5a (∆r ) 0.1 Å, ∆φ ) 1°), to produce a matrix representing the distribution of isotropic density. Using it as a simulated empirical distribution and applying the cone correction method yields a distribution as illustrated in Figure 5b. From Figure 5b, it is clear that the cone correction method only normalizes the angular effect. Figure 5a illustrates what an uncorrected distribution would look like, if no interaction whatsoever occurred. It is important to note that for Figure 5a,b, the minimum population value (violet) has a limit of zero, representing a lack of contacts at very short distances and at linear contact orientations. To calculate the elements of the isotropic corrected excess expt density distribution, Dcorr is the empirical number of ij , if nij contacts sampled during the contact count:

Dcorr ) ij

nexpt - niso ij ij niso ij

(11)

A closer look at the range of Dcorr is instructive. Even ij though it is not illustrated, applying this correction to the distribution in Figure 5a will produce a distribution yielding a figure of a single hue, violet, because the correction normalizes both angular and distance effects. However, since this correction is an excess density, the limit of the minimum values (the violet color) for this distribution is not zero. At distances very close to the hydrogen atom, no contacts are found, and therefore nexpt ) 0 in that region. However, because the ij isotropic density is always positive (it cannot be zero), it follows expt from eq 11 that Dcorr ) -1 (since niso ) 0. ij ij > 0) where nij corr Thus, it follows Dij g -1 and violet in all the isotropic density corrected rainbow contour maps (IDC RCMs) represents an excess density of -1. At distances very far from the hydrogen atom, the contact distribution becomes ideal (and therefore isotropic), and it follows that nexpt ) niso ij ij , and so corr Dij ) 0. It is not possible to analytically predict an upper limit to Dcorr for a given interaction because if Dcorr > 0 it only ij ij implies there are more contacts in the localized region, Vij ) ∆r∆φ, than expected from a purely isotropic distribution, and if Dcorr ∈ [-1,0], there are less interactions in that localized ij region than expected. An excess density is calculated to ensure

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 647

Figure 5. Rainbow contour maps of a simulation of a D-H‚‚‚A interaction: (a) a perfect noninteracting isotropic density as calculated using the isotropic density method, eq 9 and (b) applying the cone correction to the isotropic density distribution, eqs 2 and 3. it approaches zero at distances very far from the hydrogen atom, Dcorr f 0 as r f ∞. ij The IDC RCMs (with color varying as Dcorr varies) shown ij in this article were thus all obtained by first counting the contacts to obtain a distribution of contacts in a “grid” of distances and angles. Once the uncorrected distribution has been found, it is possible to do either the cone correction or the IDC (however, as shown above, the cone correction is actually implied by the IDC). Applying the IDC results in an excess density distribution and it is the excess density distribution, Dcorr ij , that is used for analysis. However, to show how the cone correction only normalizes angular effects, while the IDC normalizes contacts for both distances and angles, both correction methods were applied to actual data and will be illustrated later for some interactions. (b) Distances of Highest Incidence. To determine the distance of highest incidence, ropt, for each interaction, the is defined as distance-associated density vector Dcorr i

ERF(q) )

∑D

corr ij

j)1

N

(12)

where N is the total number of columns (thus pertaining to angles) in the IDC matrix. The vector so obtained is a representation of how excess densities concentrate at specific



2 π1/2

q

0

e-x dx 2

(13)

Once the distance-associated density is obtained, the data are subjected to a least squares approximation (from now on abbreviated LSA), optimizing the function below for the data available:

Eopt(r,r1/2,r′1/2) )

N

Dcorr ) i

distances regardless of the angles (notice that the j index is dropped since angular distribution is now averaged). These vectors can be used to determine density maxima for distance distributions numerically. The error function, ERF(q), is defined as the area under a Gaussian distribution measured from the x-axis (x ) 0) to a point q (x ) q) on the x-axis, and is defined31 mathematically as

ERF[π1/2r′1/2(r - r1/2)] - 1 2

(14)

where Eopt is the optimized error function. The form of the function (the right-hand side of eq 14) is such as to ensure Eopt ) -1 at short distances and Eopt ) 0 at longer distances in the as shown in Figure 6, and thus equivalent to Dcorr i absence of any interaction. During the LSA, the parameters (r1/2, r′1/2) are optimized. First, r1/2 describes the distance at which the distribution reaches half-height, at that distance

648

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

van den Berg and Seddon

Figure 6. The error function (eq 14) is used to account for isotropic contact distributions very close to and far away from the hydrogen atom.

Figure 7. The weighted Gaussian function (eq 15) is used to describe the accumulation of contacts where the interaction occurs. r ) r1/2 and the first derivative of the curve equals r′1/2, more formally E′opt(r ) r1/2,r1/2,r′1/2) ) r′1/2. Mathematically, the Gaussian function is defined31 as G(x) ) e-x2. For the purpose of this article, the function is redefined and termed the weighted Gaussian function:

( (

Gopt(r, h, v1/2, ropt) ) h exp ) h exp -

)

4 ln(2) (r - ropt)2 v21/2 2.773 (r - ropt) v21/2

)

2

(15)

A summary of these variables is shown in Figure 7. The constant 4 ln(2) ) 2.773 ensures that the denominator of the expression inside the exponential function equals the square of the half-width, v1/2 (thus making calculations much simpler). The variable h describes the total height of the peak and ropt the symmetry axis, or distance of highest incidence. It is hypothesized that Dcorr is equivalent to the sum of the i weighted Gaussian function and the error function, formally Dcorr ≡ (Eopt + Gopt), near the distance where the highest i number of contacts is observed. The occurrence of contacts from the second interaction shell is not modeled using this method.

(c) Details Concerning the LSA Calculations. All the LSA calculations were performed in the program KyPlot.30 There was some difficulty obtaining reproducible results (due to divergence of the parameters to be optimized in a specific LSA calculation), and a systematic method for the determination of ropt was developed. Remember, the equivalent approximated function will be the sum of the weighted Gaussian ≡ Eopt + Gopt). However, it can and error functions (Dcorr i easily happen that the weighted Gaussian function overshadows the error function (Gopt . Eopt). Due to this difficulty, it obtained must first be apwas found that the vector Dcorr i proximated using only the error function part of the expression. Doing so will ensure that the best fit for the error function for the data points in the vector will be obtained, and the weighted Gaussian peak will simply be “overlaid” onto the error function. For this method to work, it is necessary to specify good starting-point values, but that is easily accomplished by inspection. When solving the error function, it was necessary to specify an estimated starting-point value for r1/2, using a graph of Dcorr vs r: a starting point was chosen to i be close to the half-height of the error function. After choosing a suitable value for r1/2, the LSA iterative routines quickly converged to find the best values for r1/2 and r′1/2. After adding the weighted Gaussian expression to the error function expression, convergence was only obtained if an approximate starting-point for the value of ropt was given (again easily found by inspection since ropt is the symmetry axis of the first peak). After overlay of the weighted Gaussian function, the values of r1/2 and r′1/2 usually changed slightly to ensure that the final function is a LSA of all the variables. However, because the overlay of the Gaussian function often overshadowed the error function, the optimization process could not find accurate (and often not reproducible) values for r1/2 and r′1/2. Due to this restriction, these values were not used during analysis. Therefore, only the values obtained from optimization of the Gaussian function were used in the analysis. It might be possible to find better models and therefore better approximations for these distributions in future, but this will be left as a project for further research. (d) Visualization of Results. There are two intrinsic difficulties in the use of the cone correction or isotropic density correction when analyzing the data visually. First, since φ f π w sin(φ) f 0 (π ≡ 180°), division by that term in eq 2 causes asymptotic behavior near linear D-H‚‚‚A orientations. Even though the IDC does not explicitly divide by a sine term, it does so implicitly by dividing by a difference of cosines after integration, and thus the same behavior can be expected. This difficulty is sometimes also encountered in NMR research: when the concentration of the product is very low, the solvent peaks may grow large so fast that the product peaks fade into the background noise. In NMR studies, to overcome this problem, the scale of the graph is changed until the product peaks are clearly visible, while the solvent peaks are by then completely off-scale. In this article, the same approach is followed in the IDC RCMs to ensure that detail far from the asymptote is also discernible, resulting in a red-shift of color contours as shown in Figure 8a. It could be argued that such a shift in color contours will bias weaker interactions, but as shown in Figure 8b this is not the case; it only serves to clarify such effects. Since the purpose of this article is to show trends across a wide range of interactions, it was decided that a clear visual description is more informative; however, for the insistent reader a copy of the graphs before red-shifting is also available from the authors. As for NMR analysis, the choice of scale is somewhat arbitrary, since different cases will have separate rescaling requirements. The same is true for these graphs, and the upper color (red) limits to Dcorr are shown in ij Table 2. Second, the use of color contour maps neglect visual information about the actual height of the peak which becomes visible if the contour map is rotated to produce a threedimensional graph (as shown in Figure 8c). Therefore, it is difficult to compare RCMs of different interactions directly with each other. This is a serious restriction of the RCM

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 649

Table 2. The Scale Upper Limits to Dcorr for Each ij Interaction Typea interaction

upper limit

interaction

upper limit

RH2C-H‚‚‚H-CH2R O-H‚‚‚ClN-H‚‚‚ClC-H‚‚‚FC-H‚‚‚ClC-H‚‚‚BrC-H‚‚‚IC-H‚‚‚F-MZn C-H‚‚‚Cl-MZn C-H‚‚‚OdNZn C-H‚‚‚OdSZn C-H‚‚‚OdPZn C-H‚‚‚OdCZn C-H‚‚‚OdMZn

0.8 406.7 119.3 35.9 11.7 10.9 9.4 22.6 6.0 4.7 4.8 8.5 4.2 10.7

C-H‚‚‚NtCZ C-H‚‚‚OR2//OdC C-H‚‚‚SR2//OdS C-H‚‚‚F-R C-H‚‚‚Cl-R C-H‚‚‚Br-R C-H‚‚‚I-R C-H‚‚‚Br-MZn C-H‚‚‚I-MZn O-H‚‚‚OdNZn O-H‚‚‚OdSZn O-H‚‚‚OdPZn O-H‚‚‚OdCZn O-H‚‚‚OdMZn

6.8 2.5 3.0 2.2 3.0 5.3 5.1 8.5 5.3 121.6 226.6 412.2 286.3 122.1

a M ) transition metal, Z ) any substituent (the “n” indicates that the number of substituents were not fixed), R ) Z starting with a carbon atom, A//B ) all aggregates of A excluding all aggregates containing B. Note that the values were chosen arbitrarily to make visualization easier.

visualization method, but it can be solved in the graphs of the distance-associated density vectors. A comparative analysis of the height of the LSA approximated Gaussian functions for different interactions will be important, since it will be shown that there is a general trend observed between these values and the peak half-widths as obtained from the LSA calculations. Furthermore, this correlation is interesting since it may shed some light as to whether specific interaction types can be placed in a spectrum of interaction strengths.

Results and Discussion To determine whether the methods derived above could be used to qualify an interaction as a “hydrogen bonded interaction” or “van der Waals interaction”, it is necessary to apply them to interactions for which we have good experimental evidence to prove the nature of the interaction of interest. The first two interaction classes in this section have experimentally been shown to be at the extremes of the possible interactions. The Classical van der Waals Interaction. The RCMs obtained for the interaction between two methyl hydrogen atoms, RH2C-H‚‚‚H-CH2R, are illustrated in Figure 9. If the three H atoms from one methyl group were in close proximity to three H atoms from another methyl group, a total of nine possible contacts may resultsall possible contacts were counted. Exclusion of some contacts as opposed to others (e.g., only choosing the shortest contact from the nine possible contacts described earlier) would bias the counting process and make the use of the IDC method pointless. Compare Figure 9a,b with the results obtained in Figure 5a,bsthe similarities are truly remarkable. The two RCMs appear to indicate that there is absolutely no region where the contact density is higher than in the ideal case where the contacts are perfectly isotropic. The third RCM (Figure 9c) is the IDC RCM. Imagine that Figure 9a was “baselined” using Figure 5a and that gives a conceptual idea of how the correction method normalizes both angular and distance effects: if Figure 9a was perfectly isotropic, then Figure 9c would have been only one color, a violet hue that constitutes the ideal isotropic density, Dcorr ) 0. However, at roughly 2.9 Å a band of ij slightly higher density occurs, but this band is spread across all angles and therefore shows no directionality,

a feature that implies van der Waals character. Furthermore, the band shows maximum density at a distance that is much greater than the van der Waals H cutoff distance (2.9 Å > 2.40 Å ) rH w + rw ). Another characteristic of this band is that it is broadened over the region r/Å ∈ [2.5,3.5]. For the sake of continuity, it is also possible to visualize the distributions using radial plots of the IDC RCM, a visualization method that was used previously12 by this group (even though it was only applied to conecorrected data). This visualization method is illustrated in Figure 10. The region that was plotted in Figure 10a,b is shown in Figure 9c. Since the radial plots do not show more information than the RCMs, they will only be shown for a few select cases to illustrate that the results obtained in this article agree with previous research, but for most cases they will be neglected in the article and will be available as Supporting Information. Determination of ropt is shown in Figure 11. The parameters (ropt, v1/2) were found to be (2.86 Å, 1.234 Å) after LSA. Clearly, ropt was much larger than the van der Waals cutoff distance and the peak is very broad (with v1/2 greater than 1.0 Å). Classically Strong Hydrogen-Bonded Interactions. Interactions in this class are represented as Y-H‚‚‚Cl-, where Y ) O or N. It is important to compare the results obtained for van der Waals type interactions with classically strong hydrogen bonded interactions to observe the fundamental differences between the RCMs. For this purpose, the O-H‚‚‚Cl- and N-H‚‚‚Cl- interactions were analyzed. For the O-H‚‚‚ Cl- interaction, the RCM resulting from the empirical distribution is shown, Figure 12a, to illustrate its difference in comparison to the IDC RCM, Figure 12b. Also compare Figure 12 with Figure 9sthe differences are very clear. The cone-corrected RCM is not shown and will not now be further used. Figure 12 shows the RCMs of the interaction, and Figure 13 shows the related radial plots obtained from the IDC RCM. Two changes in the RCMs going from the empirical to the corrected RCM are particularly noteworthy: (a) peaking contact density shifts toward 180° (the “red-spot” shift) and (b) random isotropic noise far from the region of highest incidence fades away (random bluish blotches in the upper left-hand corner vanish). Noting that all three characteristics of van der Waals contacts are eliminated provides proof that a strong hydrogen bond is present: the distance of highest incidence (ropt ) 2.23 Å) is considerably shorter than Cl 8 the sum of the van der Waals radii (rH w + rw ) 2.96 Å), the region of highest incidence is relatively narrow (v1/2 ) 0.366 Å) and the directionality is very clear from the IDC RCM. These results are mirrored in the RCM obtained for the N-H‚‚‚Cl- interaction type shown in Figure 14 (the RCMs of the empirical distributions will not be shown for this and further interaction types). Figures 13 and 15 are the related radial plots for the O-H‚‚‚Cl- and N-H‚‚‚Cl- IDC RCMs, respectively. They mirror the observations made for the IDC RCMs. Determination of the distances of highest incidence for interactions in this class is shown in Figure 16. They are very different from the LSA graph of the van der Waals interaction type (Figure 11). Notice that this

650

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

van den Berg and Seddon

Figure 8. The effect of rescaling or rotation on the IDC RCMs; (a) the O-H‚‚‚Cl- and (b) the C-H‚‚‚Cl-R interaction rescaled and (c) the C-H‚‚‚Cl-R interaction rotated. Even though the color scale is much clearer in the second RCMs for (a) and (b), the observed trends remain the same. In (c) the masking of height information is clearly illustrated.

interaction class has very prominent Gaussian character; indeed, the error function character is almost completely hidden due to the large ratio of the Gaussian peak compared to the error function. The weighted Gaussian functions peak at 2.23 Å with a half-width of 0.366 Å for O-H‚‚‚Cl- and 2.24 Å with a half-width of 0.367 Å for N-H‚‚‚Cl-. The peaks are thus very close to each other and the peak height, h, for O-H‚‚‚Cl- is

roughly twice that of N-H‚‚‚Cl-, while the half-widths are almost equal. This interaction class provides a good example of why the parameters optimized for the error functions are difficult to reproduce. Consider again the NMR analogy, as described previously: the ratio of the maximum peak height to the background noise determines the digital observability of low-concentration products. In this

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 651

Figure 10. Radial plots with axes assigned to (a) the distance distribution in Å and (b) the angular distribution in degrees, of the RH2C-H‚‚‚H-CH2R interaction type. The color varies with the remaining parameter (angle in the distance distribution) using a rainbow color map, and Dcorr is plotted against ij the height of the axes.

Figure 9. Rainbow contour maps (RCMs) of the RH2CH‚‚‚H-CH2R interaction type: (a) the empirical, (b) the conecorrected and (c) the isotropic density corrected rainbow contour map (IDC RCM) with the classical van der Waals H cutoff distance (rH w + rw ) 2.40 Å) as the dotted line and the window (dash-dot-dot line) further described using radial plots (cf. Figure 10).

interaction class, the error function is so small compared to the Gaussian function that the information obtainable for the error function becomes meaningless (or the error in the information is too large to give reproducible values); it is completely overshadowed. Last, notice that the scale of Figure 16 is almost a 100 times that of Figure 11. This fact will be very

Figure 11. Determination of the least squares approximation (LSA) for the RH2C-H‚‚‚H-CH2R interaction class. The curve is the best fit to the data points.

important later and illustrates the inability of the IDC RCMs to convey information about the relative maxima of the IDC matrixes. Interactions with Halide Anions. The first group of weaker interactions analyzed involve the halide anions, represented as C-H‚‚‚X-, where X ) F, Cl, Br, or I. The C-H‚‚‚Cl- interaction has been analyzed previously by this group,12 and it was concluded to be a “real” hydrogen bond. This trend is confirmed as shown in Figure 17a. Radial plots are not shown since they reiterate effects already observed in the IDC RCMs and

652

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

van den Berg and Seddon

Figure 13. Radial plots of (a) the distance distribution and (b) the angular distribution of the O-H‚‚‚Cl- interaction type. (For a clear description of these radial plots, cf. Figure 11).

Figure 12. RCMs of the O-H‚‚‚Cl- interaction type: (a) the initial uncorrected RCM and (b) the IDC RCM. The IDC RCM also shows the classical van der Waals cutoff distance (rH w + rCl w ) 2.96 Å) and the window further described using radial plots (cf. Figure 13).

LSA graphs, but they are available as Supporting Information. The region of highest incidence is focused near 180° and slightly below the van der Waals cutoff distance. The difference between this IDC RCM and those for strong interactions is a “tailing” observed due to a spreading of weaker interactions away from the linear orientation in this interaction class. The database search was extended to include F, Br, and I anions. The C-H‚‚‚F- interaction (shown in Figure 17b) was particularly hard to analyze due to a lack of contacts available (118 entries) in the CSD (cf. Table 3). Even though F is accepted as a very hard atom, the IDC RCM still showed features similar to those of C-H‚‚‚Cl-; most importantly, the area of highest density occurs below the van der Waals cutoff distance! For all interactions of which the RCMs are not shown explicitly in this article, they are accepted to be similar to the RCMs shown explicitly (or previously), and are available as Supporting Information. The C-H‚‚‚Br- interaction is similar to C-H‚‚‚Cl-, but tailing becomes more pronounced as the Br anion is considerably more polarizable than its smaller cousins. This trend continues for the C-H‚‚‚I- interaction, and the region of highest

Figure 14. IDC RCM of the N-H‚‚‚Cl- interaction type showing the classical van der Waals cutoff distance Cl (rH w + rw ) 2.96 Å) and the window further described using radial plots (cf. Figure 15).

incidence dispersed away from 180 toward 160° and much closer to the classical van der Waals cutoff distance. The LSAs for this interaction type are shown in Figure 18. The F- interaction is shown separately since it is much weaker and peaks at a distance that requires a different scaling of the distance axis. The data points obtained for F- are scattered because of a lack of contact information. The remaining anion interactions are shown together and the trend observed is an increase in ropt as the size of the atom increases. The half-width for the C-H‚‚‚F- interaction is slightly larger than the rest at 0.823 Å, whereas the average half-width for the

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 653

Figure 15. Radial plots of (a) the distance distribution and (b) the angular distribution of the N-H‚‚‚Cl- interaction type.

Figure 17. The representative IDC RCMs of interactions in the C-H‚‚‚X- interaction class; (a) the C-H‚‚‚Cl- interaction and (b) the C-H‚‚‚F- interaction showing their respective van Cl H F der Waals cutoff distances (rH w + rw ) 2.96 Å, rw + rw ) 2.70 Å). From (b) the lack of contact information is clearly observable. Table 3. Number of Contacts Found for Each Interaction Type and Their Respective Compactness Constantsa interaction

Figure 16. The ropt determination LSAs for classically strong interactions. The curves represent the best fits to the datapoints of each interaction.

other three interactions is vj 1/2 ) 0.551 Å (standard deviation, σs ) 0.019). Interactions with Halocarbons. This class includes interactions that can be represented by the molecular aggregates C-H‚‚‚X-R, with X ) F, Cl, Br, or I. The C-H‚‚‚Cl-R interaction, as analyzed by Thallapally and Nangia,32 was described by them as a van der Waals interaction; recently, Nangia33 refers to this interaction type as very weakly interacting. Previously, this group12 studied this interaction and showed the inefficacy of the van der Waals cutoff distance. Recently,34 the C-H‚‚‚ Cl-R interaction was studied using Hirschfeld fingerprinting techniques as applied to a few specific interacting molecules, and the interaction was clearly observable using this method. Furthermore, hydrogen

RH2C-H‚‚‚ H-CH2R O-H‚‚‚ClN-H‚‚‚ClC-H‚‚‚FC-H‚‚‚ClC-H‚‚‚BrC-H‚‚‚IC-H‚‚‚F-MZn C-H‚‚‚Cl-MZn C-H‚‚‚OdNZn C-H‚‚‚OdSZn C-H‚‚‚OdPZn C-H‚‚‚OdCZn C-H‚‚‚OdMZn

Ntot 654971

c

interaction

0.367 C-H‚‚‚NtCZ

Ntot 109345

c 3.643

3717 108.818 C-H‚‚‚OR2// 206706 1.126 CdO 88489 2.573 8642 55.197 C-H‚‚‚SR2// OdS 872 2.198 C-H‚‚‚F-R 146182 1.811 52334 6.994 C-H‚‚‚Cl-R 182384 3.228 20174 8.195 C-H‚‚‚Br-R 41955 3.555 20734 6.395 C-H‚‚‚I-R 8486 3.212 5532 4.266 C-H‚‚‚Br-MZn 26269 5.785 175900 5.485 C-H‚‚‚I-MZn 24532 4.663 167131 3.364 O-H‚‚‚OdNZn 4741 23.169 146862 3.511 O-H‚‚‚OdSZn 9530 76.463 26659 3.736 O-H‚‚‚OdPZn 4095 135.237 681588 3.341 O-H‚‚‚OdCZn 51094 85.367 34373 4.280 O-H‚‚‚OdMZn 13231 24.162

a M ) transition metal, Z ) any substituent (the “n” indicates that the number of substituents were not fixed), R ) CZn, A//B ) all aggregates of A excluding all aggregates containing B.

bonds in this class were used to drive the complexation properties of nanoscale molecular baskets.35 Comparison of the IDC RCM obtained for this interaction with

654

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

Figure 18. The LSAs for (a) the interactions in the C-H‚‚‚ X- class excluding the C-H‚‚‚F- interaction and (b) the C-H‚‚‚ F- interaction showing a lack of continuity due to a lack of available aggregates in the CSD. Notice that the C-H‚‚‚Finteraction LSA is plotted to a slightly smaller scale on the r-axis; otherwise convergence of the LSA was difficult.

previous interactions shows a trend toward weak hydrogen bonding, as shown in Figure 19a. The directionality for these interactions is much weaker than the C-H‚‚‚X- interactions as shown in Figure 17a, and disperses down to 120°. There is also a considerable broadening of the region of highest incidence. These contacts appear at a distance slightly above the classical van der Waals contact distance, probably since the electrostatic potential is shallow on the attractive side, but they are still demonstrable interactions, not van der Waals contacts. Atoms can thus reside at longer distances in the crystal with little energy penalty. For the C-H‚‚‚F-R contact shown in Figure 19b (in contrast to the C-H‚‚‚F- contact shown in Figure 17b), enough interactions were available to do a reliable statistical analysis (3874 entries). It is similar to the C-H‚‚‚Cl-R interaction occurring closer to the arbitrary van der Waals cutoff with similar spread and broadening. The localization of contacts close to 180° is much more pronounced than for the C-H‚‚‚Cl-R interaction. This would seem to suggest a stronger directionality, but the broad tail suggests this is purely an effect of the correction method, since the sine-reciprocal tends to infinity close to 180 degrees (a singularity), the correction methods can overemphasize observed interactions close to the singularity. The IDC RCMs of the C-H‚‚‚

van den Berg and Seddon

Figure 19. The IDC RCMs representing the C-H‚‚‚X-R class; (a) the C-H‚‚‚Cl-R and (b) the C-H‚‚‚F-R interactions with their respective van der Waals cutoff distances shown Cl H F (rH w + rw ) 2.96 Å, rw + rw ) 2.70 Å).

Br-R and C-H‚‚‚I-R are similar to C-H‚‚‚Cl-R, but the C-H‚‚‚I-R interaction has a lack of contact information to be as clear as the other interactions in this class. The LSAs for these interactions are shown in Figure 20. It shows the same trend as the LSAs for the C-H‚‚‚X- interaction class with slightly broader peaks. As previously observed for the C-H‚‚‚X- interaction class, the C-H‚‚‚F-R half-width is slightly larger at 0.603 Å with the average half-width of the other three at vj 1/2 ) 0.589 Å (σs ) 0.001). Comparison of Figure 9c with Figure 19b makes quantitative interpretation clearly difficultsvisually the distinction between van der Waals contacts and very weak hydrogen bonds is difficult to resolve. To appreciate these visual similarities, the weak hydrogen-bonded interactions must be interpreted as a combination of electrostatic and dispersive behaviors. However, comparison of Figure 11 with Figure 20a is noteworthy, and illustrated, for the C-H‚‚‚F-R interaction, in Figure 20b. Specifically, the C-H‚‚‚F-R peak, with half-width of roughly 0.6 Å, is much narrower than the RH2CH‚‚‚H-CH2R peak with half-width of roughly 1.2 Å. Later, the information obtainable from the LSA graphs will be used to analyze the variable nature of different interaction classes in greater detail (cf. Figure 28).

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 655

Figure 21. The IDC RCM of the C-H‚‚‚Cl-MZn interaction Cl (rH w + rw ) 2.96 Å) representing the general trend of RCMs in this class.

Figure 20. The LSAs for (a) interactions in the C-H‚‚‚X-R class; (b) interactions in the C-H‚‚‚F-R and RH2C-H‚‚‚HCH2R interaction class as comparison, illustrating the difference between a very weak hydrogen bond (green) and a van der Waals contacts (blue).

However, it is reiterated at this point that both the IDC RCMs and the LSA graphs have to be analyzed in tandem to get a clear picture of the localization of contacts; the IDC RCMs lack information of relative peak height, while the LSA graphs lack directional information. Interactions with Halides Bound to Transition Metals in Uncharged Complexes. Interactions in this class are represented as C-H‚‚‚X-MZn, where X ) F, Cl, Br, or I; M is any transition metal and Z any substituent (the n below Z implies that the number of substituents were not fixed). The C-H‚‚‚Cl-MZn interaction was analyzed by Thallapally and Nangia,32 and they observed that this interaction is a weak hydrogen bond. This is confirmed by the RCMs in Figure 21, but it is noteworthy that this interaction, while more directional than the C-H‚‚‚Cl-R interaction, is just as broad. This implies that the difference between the C-H‚‚‚Cl-MZn and C-H‚‚‚Cl-R interaction might not be as large as previously thoughtscertainly not big enough to arbitrarily fall into two different interaction classes (as suggested by Thallapally and Nangia).32 The distance of highest incidence for the C-H‚‚‚Cl-MZn interaction is slightly shorter in comparison with the C-H‚‚‚Cl-R interaction, at ropt ) 2.99 Å (for C-H‚‚‚Cl-

R, ropt ) 3.17 Å). The C-H‚‚‚F-MZn interaction was analyzed by Brammer et al.,36 and they showed this interaction to be hydrogen-bonded; even though it again suffers from a lack of available data in the CSD, its region of highest density is localized and occurs below the classical cutoff distance so it is still a demonstrable interaction. The C-H‚‚‚Br-MZn and C-H‚‚‚I-MZn interactions follow a similar trend but are increasingly less directional. Figure 22 show the LSAs for these interactions. As expected, they show a trend similar to LSA graphs for the other halide interaction classes. However, the half-widths in this class are all very similar and the average half-width is vj 1/2 ) 0.557 Å (σs ) 0.011). Interactions with π-Bonded Oxygen. Interactions in this class can be represented as X-H‚‚‚OdYZn, where X ) O or C; and Y ) C, N, P, S, or M. This group of interactions has not yet been studied as a whole, using statistical methods. Due to the electron distribution around the oxygen atom, a hydrogen-bonded interaction is expected to occur.14-17,37 Every C-H interaction was compared to a similar O-H interaction. The IDC RCMs of the O-H contact interactions suggest the occurrence of strong hydrogen bonds as expected. These interactions occur well below the classical van der Waals cutoff and are directional and narrow in region. An example of these interactions is shown in Figure 23asit is representative of all contacts in the O-H contact subclass. When the IDC RCMs of the C-H contact interactions are analyzed, a high degree of directionality is still observed. The region of highest incidence is much broader and dispersed away from 180°, but altogether the graphical topography suggests stronger interactions when compared to interactions with halogens. Figure 23b is a representative example of contacts in this subclass. Figure 24 illustrates the LSA graphs for these interactions. Notice that LSAs for a variety of environments in this class remain very close to each other; this implies that the secondary atom does not appreciably influence ropt. The C-H interaction LSA graphs appear to fall almost on top of each other, and this point to interactions of similar strength. However, O-H interactions differ markedly from each other as far as peak

656

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

Figure 22. The LSAs for interactions in the class C-H‚‚‚XMZn; (a) shows only the C-H‚‚‚F-MZn interaction, whereas (b) shows the other interactions.

height is concerned. Does this point to markedly different strengths? It would appear that it is not necessarily only the peak height that is correlated directly to strength-of-interaction, but possibly also the peak halfwidth. The half-widths for both subclasses are of the same order within each subclass and the average halfwidths are vj 1/2 ) 0.371 Å (σs ) 0.028) for the O-H‚‚‚ OdYZn subclass and vj 1/2 ) 0.593 Å (σs ) 0.007) for the C-H‚‚‚OdYZn subclass; this shows that the O-H class has a much narrower spread (and therefore appears stronger) than the C-H subclass. Interactions with Ethers and Thioethers. Interactions in this class can be represented as C-H‚‚‚YR2 where Y ) O or S. Bond and Jones38 studied the C-H‚‚‚SCM interaction using the cone correction and standard scatter plot techniques. They concluded that the interaction does not have directing properties in the crystalline state. Studying this interaction class was hampered by the presence of CdO and OdS functionalities in the same aggregates (in esters and sulfoxides, for example), since hydrogen bonds were more likely to exist as C-H‚‚‚OdC than as C-H‚‚‚OR2. However, it was possible to set up the search so specifically exclude aggregates containing these functionalities. The LSA determination and IDC RCMs of these interactions are shown in Figure 25. From these graphs, it appears that the two interactions are of roughly similar strength. The

van den Berg and Seddon

Figure 23. The IDC RCMs of interactions (a) O-H‚‚‚OdSZn O and (b) C-H‚‚‚OdSZn including cutoff distances (rH w + rw ) 2.70 Å).

C-H‚‚‚OR2 interaction peaks at 2.78 Å, whereas the C-H‚‚‚SR2 interaction peaks at 3.21 Å and the halfwidths are, respectively, 0.372 and 0.524 Å. Using the IDC RCM method it was thus possible to show that the C-H‚‚‚SR2 interaction still have directional properties, by extracting the interacting contacts from random background contacts. Interactions Involving Nitrogen. Currently, only the interactions represented as C-H‚‚‚NtCZ have been studied. At first glance, it can be expected that this situation interact due to available electron density on the cyanide nitrogen atom. This is confirmed in the IDC RCM shown in Figure 26 and the LSA determination in Figure 27. The interaction has all the features expected for a hydrogen-bonded interaction; specifically, the LSA resulted in ropt ) 2.84 Å and v1/2 ) 0.644 Å. Comparative Analysis of Interactions. So far, each interaction class was analyzed separately to show trends between different specific interactions. Now all interactions will be studied as a whole. This analysis is somewhat harder as explained earlier. It is important that there should be enough data in the database to ensure a minimization of statistical error. As shown in Table 3, this is generally the case. Table 3 is a summary of the total number of contacts counted, Ntot, for each

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 657

Figure 24. The LSA graphs for (a) the O-H‚‚‚OdYZn subclass and (b) the C-H‚‚‚OdYZn subclass. Notice that the peaks in (a) are much narrower than those in (b).

interaction analyzed. Looking at the IDC RCMs in general, an interesting trend can be seen: stronger interactions tend to have highly localized regions of highest density. To quantify this localization of contacts, a compactness ratio χ, is defined in eq 16.

χ)

h v1/2

(16)

To rationalize this choice, refer to the weighted Gaussian curve shown in Figure 7. It is an approximation of the excess density of contacts at a certain distance, and therefore the concentration of contacts at the peak is directly proportional to the height, h, of the peak (since it is a measure of the excess density) and inversely proportional to the half-width (since it is a measure of the localization). The ratio χ will have larger values for either higher or narrower peaks. The values of this ratio for each case are shown in Table 3. A general, but not quantitative, observation is larger values of χ for interactions that are conceptually stronger (e.g., comparing O-H‚‚‚Cl- to C-H‚‚‚Cl-). Even though it is not a quantitative measure of strength between various confirmed hydrogen bonds, χ was less than unity only for the van der Waals contact case RH2C-H‚‚‚H-CH2R. Further study may show this to be true for all possible searches that represent pure van der Waals contacts! For all the interactions, ropt was determined and compared with classical van der Waals cutoff distances.

Figure 25. The LSA and IDC RCMs for interactions in the class C-H‚‚‚YR2; (a) the LSA determination and RCMs for (b) O C-H‚‚‚OR2 and (c) C-H‚‚‚SR2 interactions; (rH w + rw ) 2.70 Å, H S rw + rw ) 3.03 Å).

Since χ is a measure of the compactness of interactions near the distance of highest incidence, it was hypothesized that there might be a correlation between χ and ropt. To highlight weaker interactions, it is possible to use a log-axis for plotting χ-values vs ropt (such an axis system generally accentuates differences in small numbers). Figure 28 is an overview of the results obtained: notice the segregation of interactions into various regions that correspond to varied strengths. This can

658

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

Figure 26. The IDC RCM of the C-H‚‚‚NtCZ interaction N (rH w + rw ) 2.80 Å).

van den Berg and Seddon

Figure 28. The relationship between ropt and log10 (χ). The log-axis was chosen to accentuate differences is the small χ-values obtained for weak interactions. van der Waals contacts are shown as diamonds ([) and hydrogen-bonded interactions as circles (b). Regions of different conceptual interaction strength were highlighted using different background colors: peach-red indicates a strong interaction, yellow indicates a weak interaction, green indicates a very weak interaction, and magenta indicates a van der Waals contact. Circles with similar colors represent specific contacts in the same class as shown by the labels connected to the circles, e.g., all circles filled with a cyan color fall in the class Y-H‚‚‚Cl- (cf. Tables 3 and 5). The scale on the right depicts interaction energies calculated using gaseous phase super-molecule methods, previously assigned by various authors.4,39-41 This scale does not reflect on the actual interaction energies in the crystalline state, but is included to illustrate how calculated gas-phase energies roughly correlate to different χ-values.

Figure 27. The LSA determination of the C-H‚‚‚NtCZ interaction.

be seen, in a simplified form, in Figure 29, where the data points have been excised. An interesting observation is the changes of ropt for the halogen-containing interaction classes C-H‚‚‚ X-, C-H‚‚‚X-MZn and C-H‚‚‚X-R. Define the value ∆ropt(A f B) ) ropt(B) - ropt(A), where A and B are two different interaction classes (e.g., C-H‚‚‚Cl- and CH‚‚‚X-MZn). These values are illustrated in Table 4. It is noteworthy that these changes are of the same order within each group (as shown in the columns of the table) and mirrors the weakening of these interactions, as less electron density becomes available for an interaction to occur. Since evaluating the validity of the van der Waals cutoff distance was the purpose of this analysis, distances of highest incidence were compared to those proposed by Bondi.8 In Table 5, these results are listed. Notice that Bondi’s assumption that the distances of highest incidence are equal for similar interactions only differing in atomic environments is clearly violated as shown in Figure 30a. For very strong interactions, ropt lies well below the cutoff distance, but for weaker interactions, the separation is not as clear. For some weak interactions, distances of highest incidence are observed slightly above the artificially proposed cutoff.

Figure 29. A simplified illustration of Figure 28, illustrating the segregation of contacts into different interaction strengths by removing the data points of each interaction. As shown in Figure 28, regions of different conceptual interaction strength were highlighted using different background colors: peachred indicates a strong interaction, yellow indicates a weak interaction, green indicates a very weak interaction, and magenta indicates a van der Waals contact. Table 4. Changes in the Distances of Highest Incidence for the Halogen-Containing Interaction Classes halogen, X

∆ropt(X- f X-MZn)/Å

∆ropt(X- f X-R)/Å

F Cl Br I

0.075 0.070 0.098 0.082

0.383 0.253 0.250 0.221

Thus, Steiner and Desiraju’s proposal is verified; the variability of intermolecular effects due to the specific atomic environments cannot be oversimplified by the

C-H‚‚‚X Hydrogen Bonding in Crystalline State

Crystal Growth & Design, Vol. 3, No. 5, 2003 659

Table 5. A Comparison between ropt for Each Interaction and the Proposed Cutoffsa interaction O-H‚‚‚Cl-

N-H‚‚‚ClC-H‚‚‚FC-H‚‚‚F-MZn C-H‚‚‚F-R C-H‚‚‚ClC-H‚‚‚Cl-MZn C-H‚‚‚Cl-R C-H‚‚‚BrC-H‚‚‚Br-MZn C-H‚‚‚Br-R C-H‚‚‚IC-H‚‚‚I-MZn C-H‚‚‚I-R

ropt/Å

Z rH w + rw/Å

interaction

ropt/Å

Z rH w + rw/Å

2.231 2.238 2.420 2.494 2.803 2.916 2.987 3.169 3.025 3.123 3.276 3.269 3.351 3.490

2.96 2.96 2.70 2.70 2.70 2.96 2.96 2.96 3.05 3.05 3.05 3.16 3.16 3.16

O-H‚‚‚OdNZn O-H‚‚‚OdSZn O-H‚‚‚OdPZn O-H‚‚‚OdCZn O-H‚‚‚O)MZn C-H‚‚‚OdNZn C-H‚‚‚OdSZn C-H‚‚‚OdPZn C-H‚‚‚OdCZn C-H‚‚‚O)MZn C-H‚‚‚OR2//CdO C-H‚‚‚SR2//OdS C-H‚‚‚NtCZ RH2C-H‚‚‚H-CH2R

1.947 1.882 1.783 1.821 1.886 2.715 2.664 2.594 2.680 2.657 2.775 3.211 2.838 2.860

2.70 2.70 2.70 2.70 2.70 2.70 2.70 2.70 2.70 2.70 2.72 3.03 2.80 2.40

a Italicized contacts are van der Waals interactions. Z refers to the second atom needed to calculate the van der Waals cutoff distance, cf. Table 1.

(shown as open dots) are compared to the distances of highest incidence determined during this analysis. In Figure 30b, the H‚‚‚O interaction, as an example, includes all interaction types in which those two atoms were involved. The stronger interactions occurred below the cutoff distance, but weaker interactions also occurred slightly above that distance. A general effect observed in all the LSAs for cases that were hydrogen bonded is that beyond the peak modeled by the Gaussian curve, a region is found where there are less contacts than the isotropic density. It appears the existence of an interaction generally causes contacts to concentrate in a certain region, and to achieve this end, lowers the chance of finding contacts nearby that region. Currently, this “sucking” effect has not been modeled explicitly, but it is noteworthy that this effect was only observed for hydrogen-bonded interactions, first because it gives another general proof-of-existence using only the LSA graphs. Second, this effect may give an insight into the strength of the interaction if it could be modeled. Such a model is left for future research, but together with Gaussian peak height and half-width, there is now yet another feature of these graphs that may become important for statistical analysis. Conclusion

Figure 30. Comparison of van der Waals radii and distances of highest incidence calculated earlier; in (a), the van der Waals cutoff distance is shown as a solid line at 2.96 Å, while ropt for each interaction is shown as the open dot and the half-width of the peak is shown in the error bar; in (b) a spread of distances of highest incidence for all interactions is shown, the open dots show the distances determined using Bondi’s8 van der Waals radii, whereas the bar shows the distribution of the distances of highest incidence obtained for each interaction type in this article over all analyzed atomic environments.

use of the previously proposed van der Waals cutoffs. Figure 30b is perhaps a clearer reflection of Table 5; the various van der Waals radii described by Bondi8

Many authors, including Pauling,1 have developed schemes to define the van der Waals radii of atoms. Bondi’s8 method was pragmatic, and his results were reproducible, leading to a verifiable set of values enabling researchers to calculate molecular volumes (i.e., the volume impenetrable by other molecules). However, when this method is used to establish the occurrence of hydrogen bonds, there is a tendency to forget why the van der Waals radii were originally determined! The van der Waals radius is supposed to assist a researcher in density calculations (through which the van der Waals radii were verified) and for this purpose, the van der Waals radii are very useful indeed. However, this investigation shows that hydrogen bonding phenomena occur across a wide range of C-H interaction classes and should not be subjected to the “van der Waals radius test” by which an arbitrary cutoff value determines the hydrogen bonding character of the interaction involved. Using a statistical method, it was possible to show how a myriad of separate molecular aggregates orient

660 Crystal Growth & Design, Vol. 3, No. 5, 2003

hydrogen-containing functionalities in a fashion that approaches a linear orientation. Because of the forces at play during the process of crystal packing, this preference for a linear orientation is not immediately evident, but the development of the isotropic density correction method made this easier to observe by removing the occurrence of random isotropic contact density from the data. The directional property of all the contacts analyzed is the key factor necessary to distinguish hydrogen-bonded interactions from van der Waals contacts. After showing that most of the interactions analyzed do in fact have linear preference, distances of highest incidence were calculated based on simplistic, but still effective, modelssusing a least squares approximation method. By application of this method, all interactions that showed linear preference could be compared to the sum of van der Waals radii and even though some stronger interactions did have highest contact densities below the cutoff distance, this was not always the case, especially for very weak interactions. Using the compactness ratio, it was possible to place the different contact classes in a qualitative spectrum of interaction strengths. Do not let this confuse though; this technique is valid only when looking at interactions in general, and specific molecular aggregates do not necessarily follow this rule. To determine whether a specific contact is due to a hydrogen bond and if there is a correlation between charge distribution and the observed crystal structures, ab initio methods have to be used. Currently, there are approaches available to study differences between van der Waals contacts and hydrogen bonds using ab initio methods; Popelier and Smith recently reviewed the well-known AIM theoretical method based on a quantum topological ansatz.42 In an effort to ease the computation time of high-level computational methods, Espinosa-Garcı´a applied integrated ab initio molecular orbital-molecular orbital (IMOMO) methods to establish their applicability to variational transition-state theory studies of hydrogen bonding energies.43-46 However, since these bonds may play a very important role in the development of crystal engineering methods and their respective strengths underlie various hierarchical approaches to predict self-assembly, it is very important that we carefully establish empirically how they affect overall stability of crystalline structures. Once the occurrence of these interactions have been generalized, our predictive methods will improve considerablysbecause the hydrogen bond plays the most important role in the area of “noncovalent” chemistry. The use of the isotropic density correction method enables us to focus beyond random atomic noise that plays no role in an interaction and therefore might help further development of empirical approaches based on statistical studies. The method employed to determine distances of highest incidence is very robust and a better model would probably yield results that ultimately contain physicochemical meaning. If a model that solves both distance and angle distributions could be developed, it might be possible to incorporate such a model into modern computational methods; there might even be an overall relationship between different interaction classes,

van den Berg and Seddon

which might become apparent as soon as we find an efficient model. However, as our understanding of hydrogen bonds develops, it is clear that paradigm of the “black-and-white model” is not efficient and the philosophy of hydrogen bonding interactions has to expand to deal with a spectrum of interaction strengths depending on the electronic properties of the molecular aggregates involved. Though the controversy behind very weak hydrogen-bonded interactions may sometimes appear nonsensical, there is an important thing to remember: if the number of weak interactions in a crystalline solid is very large, and they all contribute just a small amount to the stabilization energy, then their cumulative contribution cannot be ignored; indeed, its broader influence on structure has been shown in this paper. Acknowledgment. We wish to thank Drs. Timothy A. Evans and Christer Aakero¨y for their pioneering work with the original statistical searching methods, Christof Hanke for his comments on the development of the isotropic density correction method, Dr. Stuart James for his comments on the readability of the text, and Prof. Mike Green and Dr. Alta Ranwell from the Sasol Technology Research and Development Division for their sponsorship (JAvdB, Grant R8112QLL) and support. We would also like to express our gratitude to the referees, whose invaluable comments significantly strengthened the manuscript. Supporting Information Available: A complete set of bmp images of the isotropic density corrected rainbow contour maps for all the different interaction types analyzed, as well as radial plots, in a single zip file. This material is available free of charge via the Internet at http://pubs.acs.org.

Symbols and Nomenclature CSD IDC RCM LSA rw (r, φ) lcirc

nij ncc ij W

Vinit Vcorr Vfin nfin ij (x, y, z) (r, φ, φ) Jxyz(r, φ, φ) Fjiso

Cambridge Structural Database isotropic density corrected rainbow contour map least squares approximation van der Waals radius as defined by Bondi the distance and angle parameters between a hydrogen and donor atom D-H‚‚‚A as determined by a search in the CSD circumference of the circle formed by the perimeter of the base of the cone emanating form the hydrogen atom as used in the cone correction method an element of the distance vs. angle distribution matrix, i is the distance index and j the angle index an element of the distribution matrix obtained after the cone correction method was applied, before normalization a weight factor required to ensure the volume of the cone-corrected distribution matrix is equal to the initial distribution, the normalization factor volume of the initial distribution volume of the cone-corrected distribution volume of the final distribution, equal to the initial volume after the weight factor was applied element of the normalized cone-corrected distribution coordinates in Euclidean space coordinates in spherical space the Jacobian matrix the average isotropic density

C-H‚‚‚X Hydrogen Bonding in Crystalline State Ntot Vtot (rmin, rmax) niso ij Vij nexpt ij Dcorr ij

Dcorr i ERF(x) Eopt E1/2, E′1/2 Gopt (h,v1/2,ropt) σs χ Å

the total number of sampled contacts the total sampling volume the minimum and maximum radius used to describe the sampling volume an element of is isotropic contact distribution a volume element corresponding to a single element in the initial distribution the actual number of contacts sampled from the CSD in a single volume element an isotropic density corrected excess density element, obtained after applying the isotropic density correction to the initial distribution matrix distance-associated density vector the statistic error function the optimized error function after application of the LSA the distance and slope at the half-height of the optimized error function the optimized Gaussian function after application of the LSA the total height of the Gaussian function, the half-width and the distance of highest incidence the standard deviation of a sampled data series the compactness ratio unit of measurement equals 10-1 nm.

References (1) Pauling, L. The Nature of the Chemical Bond, 3rd ed.; Cornell University Press: New York, 1960. (2) Pimentel, G. C.; McClellan, A. L. The Hydrogen Bond; Freeman and Company: San Francisco, 1960. (3) Vinogradov, S. N. Hydrogen Bonding; Van Nostrand Reinhold Company: New York, 1971. (4) Steiner, T. Angew. Chem., Int. Ed. 2002, 41, 48-76. (5) Steiner, T.; Desiraju, G. R. Chem. Commun. 1998, 891892. (6) Cotton, F. A.; Daniels, L. M.; Jordan, G. T., IV; Murillo, C. A. Chem. Commun. 1997, 1673-1674. (7) Stone, A. J. The Theory of Intermolecular Forces; Clarendon Press: Oxford, 1996; Vol. 32. (8) Bondi, A. J. Chem. Phys. 1964, 68, 441-451. (9) Watson, J.; Berry, A. DNA: The Secret of Life; Knopf, 2003. (10) Watson, J.; Gilman, M.; Witkowski, J.; Zoller, M. Recombinant DNA, 2nd ed.; W. H. Freeman and Company: San Francisco, 1992. (11) Raymo, F. M.; Bartberger, M. D.; Houk, K. N.; Stoddart, J. F. J. Am. Chem. Soc. 2001, 123, 9264-9267. (12) Aakero¨y, C. B.; Evans, T. A.; Seddon, K. R.; Pa´linko´, I. New J. Chem. 1999, 23, 145-152. (13) Moran, S.; Ren, R. X.-F. J. Am. Chem. Soc. 1997, 119, 20562057. (14) Desiraju, G. R. Acc. Chem. Res. 1991, 24, 290-296. (15) Steiner, T. Chem. Commun. 1997, 727-734.

Crystal Growth & Design, Vol. 3, No. 5, 2003 661 (16) Taylor, R.; Kennard, O. J. Am. Chem. Soc. 1982, 104, 50635070. (17) Steiner, T. J. Chem. Soc., Chem. Commun. 1994, 23412342. (18) Price, S. L.; Stone, A. J.; Lucas, J.; Rowland, R. S.; Thornley, A. E. J. Am. Chem. Soc. 1994, 116, 4910-4918. (19) Potter, B. S.; Palmer, R. A.; Withnall, R.; Chowdhry, B. Z.; Price, S. L. J. Mol. Struct. 1999, 485-486, 349-361. (20) Price, S. L.; Wibley, K. S. J. Phys. Chem. A 1997, 101, 21982206. (21) Umeyama, H.; Morokuma, K. J. Am. Chem. Soc. 1976, 99, 1316-1341. (22) Allen, F. H. Acta Crystallogr. Sect. B 2002, 58, 380-388. (23) Allen, F. H.; Motherwell, W. D. S. Acta Crystallogr. Sect. B 2002, 58, 407-422. (24) Aullo´n, G.; Bellamy, D.; Brammer, L.; Bruton, E. A.; Orpen, A. G. Chem. Commun. 1998, 653-654. (25) Aakero¨y, C. B.; Seddon, K. R. Chem. Soc. Rev. 1993, 397407. (26) Desiraju, G. R. Acc. Chem. Res. 2002, 35, 565-573. (27) Desiraju, G. R.; Ciunik, Z. Chem. Commun. 2001, 703-704. (28) Kroon, J.; Kanters, J. A. Nature 1974, 248, 667-669. (29) Bruno, I. J.; Cole, J. C.; Edgington, P. R.; Kessler, M.; Macrae, C. F.; McCabe, P.; Pearson, J.; Taylor, R. Acta Crystallogr. Sect. B 2002, 58, 389-397. (30) Yoshioka, K. KyPlot version 2.0 beta 14 (32 bit). The latest release, version 3.0, is now available for purchase from Kyence Inc. on the Internet at http://www.kyence.com. (31) Stocker, H.; Harris, J. W. Handbook of Mathematics and Computational Science; Springer-Verlag: New York, 1998. (32) Thallapally, P. K.; Nangia, A. CrystEngComm 2001, 3, 1-6. (33) Nangia, A. CrystEngComm 2002, 4, 93-101. (34) Spackman, M. A.; McKinnon, J. J. CrystEngComm 2002, 4, 378-392. (35) Gibb, C. L. D.; Stevens, E. D.; Gibb, B. C. J. Am. Chem. Soc. 2001, 123, 5849-5850. (36) Brammer, L.; Bruton, E. A.; Sherwood, P. New J. Chem. 1999, 23, 965-968. (37) Novoa, J. J.; Lafuente, P.; Mota, F. Chem. Phys. Lett. 1998, 290, 519-525. (38) Bond, A. D.; Jones, W. J. Chem. Soc., Dalton Trans. 2001, 3045-3051. (39) Rovira, M. C.; Novoa, J. J. Chem. Phys. Lett. 1997, 279, 140150. (40) Novoa, J. J.; Whangbo, M.-H.; Williams, J. M. Chem. Phys. Lett. 1991, 180, 241-248. (41) Rovira, M. C.; Novoa, J. J.; Whangbo, M.-H.; Williams, J. M. Chem. Phys. 1995, 200, 319-335. (42) Popelier, P. L. A.; Smith, P. J. In Chemical Modelling: Applications and Theory; Hinchliffe, A., Ed.; The Royal Society of Chemistry: Cambridge, 2001; Vol. 2, pp 391448. (43) Espinosa-Garcı´a, J. J. Phys. Chem. A 2000, 104, 7537-7544. (44) Espinosa-Garcı´a, J. J. Phys. Chem. A 2002, 106, 5686-5696. (45) Espinosa-Garcı´a, J. J. Phys. Chem. A 2003, 107, 1618-1626. (46) Espinosa-Garcı´a, J.; Do´be´, S. J. Phys. Chem. A 1999, 103, 6387-6393.

CG034083H