- Email: [email protected]
Molecular mechanics and molecular shape
Journal of Molecular Structure (Tbochem), 179 (1988) 83-98 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
MOLECULAR
MECHANICS AND MOLECULAR
83
SHAPE
Part VI. The response of simple molecules to bimolecular association****
A.Y. MEYER Department of Organic Chemistry, Hebrew University, Jerusalem (Israel) (Received 20 July 1987)
ABSTRACT A molecular-mechanical force field, originally developed to calculate the structure and conformation of single molecules, has been applied to pairs of identical molecules lying at close proximity. The prototypes examined are very simple - short alkanes and short haloalkanes - but do include cases with conformational diversity and possible electrostatic effects. The optimized mutual orientations include the parallel alignment of linear alkanes, head-to-tail orientations of haloalkanes, as well as the head-to-head association, known from the crystallographic literature. The effect of aggregation on the conformation of the separate partners has been studied, and the intramolecular response characterized: in all cases examined, the partners contract. The computed energy gains are of the expected order of magnitude, the main contribution coming always from what molecular-mechanists label “non-bonded attraction”. Graphical representations of optimized orientations show the “protrusion-in-cleft” matching of molecular surfaces.
INTRODUCTION
The author has reported on the molecular mechanics [3] and on the size and shape [1,2,4] of organic halogen compounds. Two trends in the current literature make it timely to go beyond single molecules, and apply the existing computational techniques to bimolecular aggregates. One is the rising interest in intramolecular geometrical changes that accompany molecular association [5]. The other is the accumulation of data regarding angular preferences of intermolecular forces [ 61. The aim here is not to re-optimize force field parameters and reproduce documented data. On the contrary, it is to examine to what extent a single-molecule force field can, together with appropriate descriptors of size and shape, contend with association and add to the conceptual background. *Dedicated to Professor Bernard Pullman. **For Parts IV and V, see refs. 1 and 2.
0166-1280/88/$03.50
0 1988 Elsevier Science Publishers B.V.
Thus, the questions we ask are primarily of principle, in the spirit of former molecular-mechanical examinations of the hexane crystal [ 7,8] and of hexaneassociates [ 91. To get answers, however, extensive numerical checking is unavoidable. As a start, pairs of very simple molecules were examined, alkanes and haloalkanes. The force field has been described [3]. Use of a different computational technique, molecular-mechanical or another, would undoubtedly affect numerical detail. It should not affect the lines of approach, nor, with evident reservations, the nature of results. The main difference of our force field from others is that all bonds (including C-H) are attributed bond dipoles, except C-C bonds when neither carbon atom carries a heteroatomic substituent. The descriptor-evaluating program has also been described [ 1,2,4]. The literature has accorded little attention to dimers of the types to be discussed. We therefore find it helpful to provide computer-drawn cuts through the van der Waals’ bodies of the dimers. To obtain these, each atom is represented by a sphere of the appropriate van der Waals’ radius [lo], each molecule thus becoming a system of interlocking spheres. To help in viewing, nuclei within 0.04 nm of the plane of cut are marked by circles, of radius equal to half the atomic radius. Nuclei farther away are not shown. Sparse hatching indicates an heteroatom. To convey a feeling for aggregate size, all figures include the “equivalent ellipse”. This is the generating ellipse of an ellipsoid of revolution that (a) exactly circumscribes the atom farthest away from the center, and (b) whose volume is equal to the volume of the van der Waals body. In most figures, the center is center of mass, and the plane shown is the first inertial plane (i.e., defined by 1, and 1s [ 111, here denoted xy). In a few instances we preferred the center of volume and the first plane by volume. These are obtained exactly as in the inertial case, except that atoms are weighted by cubed van der Waals radii. EXTREMA ON THE POTENTIAL SURFACE
As is well known, in “optimizing a chemical structure”, the investigator starts by guessing a “trial geometry”. Next, the computer program is allowed to modify the geometry by steps, such that energies computed for sequential structures go steadily down. The process stops at a “final geometry”, whereat further lowering is not feasible. Since minima are separated by ridges and the process works downhill, a minimum risks being overlooked if the trial geometry is not within the proper valley. In dealing with internal rotation in single molecules, constraints of chemical bonding, previous knowledge, and symmetry considerations [ 121 serve as guides. When in doubt, one may “drive” along the saddles [ 13,141. By contrast, little is known about preferred relative orientations of two (or more) molecules, and detailed scans of potential surfaces - allowing relaxation in all but a few degrees of freedom - are not worth the while. Consider the constellations [ 151 of Me&X-dimers (X = halogen). Me&X
is conformationally uniform and highly symmetric, and has electrostatically defined head (X) and tail (Me&). One expects a minimum when the C-X moieties of the two molecules are antiparallel to each other. Computation reveals splitting in two separate minima. One is “singly flanking”, A,: the Clatom in one sole partner is flanked by two methyl groups of the other partner (Fig. 1). The other is “doubly flanking”, AZ: the Cl-atom in both partners is thus flanked (Fig. 2). AZ is a deeper minimum and, for all we know, the minimum minimorum. Computed energy differences, E (A, ) -E (A,), are 0.8 kcal (molofpairs)-’ [16] forX=F,1.0forX=C1,Br,and0.3forX=I.A1andA, are separated by a barrier to external rotation. When a residue of local symmetry C,, rotates internally (e.g., internal rotation in Me-Me), the periodicity is 3, there are 3 minima in a cycle and 6 extrema in all. In the Fourier series V(m) = ( V,/2) (1 -cos 3~) + ( V,/ 2) (l-cos 6w)+ .**, the term in V, has the correct symmetry and is leading [ 171. By contrast, the external rotation A,*A,* ** has 6 minima in a cycle and 12 extrema in all. The term in V, does not reflect correctly the sequence of events, and calculations suggest that it is not leading. More specifically, assume that maxima occur halfway between A, and A2 (rotated through 46 with respect to either), and denote by dA and dB, respectively, the energy of Al and of the barrier above AZ. In the Fourier series, V, =dA and V, =dB - (dA/2). To estimate dB in Me&I, one of the molecules in constellation A, was rotated through 46 about C-I. The energy rose to Y
Y
L
-----
UNIT: 0.1
NM
XL_--__X UNIT: 0.1
NM
Fig. 1. “Singly-flanking” constellation (A,) of (Me&I),. In all figures the molecule is cut across the first inertial plane, except in Figs. 5 and 9 where this is the first plane by volume. All figures include the generating ellipse of the equivalent ellipsoid. Fig. 2. “Doubly-flanking”constellation
(A,) of (Me&I),.
86
2.57 kcal mol-l. Holding in place the central carbon atom and one other carbon in both units, the assembly was allowed to relax. dB fell to 0.98 kcal mol-l. If the higher estimate is the more realistic [ 81, V, - 0.3 and V, - 2. At any rate, V, is very much greater than V,. The geometry at the barrier is shown in Fig. 3. In all three constellations (Figs. l-3), the van der Waals’ body of a Me&I-unit has virtually the same size: its volume is N 0.109 nm3 molecule-‘, its surface area is N 1.35 nm2 molecule-i. What the constellation does affect are the dimensions of the intermolecular channel. This is the “bimolecular variant” of a conclusion drawn for single molecules: configurational isomers differ in density mainly because of differences in packing, less so because of differences in size [ 41. At the barrier (Fig. 3), the partners touch or interpenetrate, and the channel is blocked. It should be stressed that the energy at the barrier is still well below that of two non-interacting molecules. At least in this case, external rotation within the aggregate is cheaper to activate than dissociation followed by reorientation. Generalizing now, it may happen that normal contacts along extensive ranges of the interacting surfaces require a close contact at some point. Then, it may again happen that stabilization due to the former outweighs the destabilization caused by the latter. Now, the van der Waals’ volume of a molecule was defined in 1964 as the volume that is impenetrable to other molecules with thermal energies at ordinary temperatures [ 181. Later work does not substantiate this view. Recent scans of the Cambridge Crystallographic Data Base revealed 729 43 =
l______x UNIT:
0.1
NM
L-_-_-X UNIT:
0.1
Fig. 3. Barrier separating constellations A, and A2 of (Me,CI),. Fig. 4. “Head-to-head” constellation
(HH) of (Me&Cl),.
NM
instances of interpenetrating C-X. * (0,N) contacts [6] and 86 instances of interpenetrating C-1. **(O,N,S) contacts [ 191. An intramolecular instance is provided by 2,2 ’ -dihalobiphenyls, where halogen atoms approach closer than the sum of theirvan der Waals’s radii [ 201. Our simple prototype shows that intdr!enetration does not imply an energy-rise. It may well enable a considerable energy-gain. Two other minima were located on the potential surface of (Me&X)2, this time separated from A, (Fig. 2) by barriers to translation. In one, the partners are head-to-tail (HT, Me&X* - *Me&X); in the other, they are head-to-head (HH, Me&X* * *XCMe,). The more interesting constellation, HH, is shown in Fig. 4. By computation, it lies 1.3 kcal mol-’ above A2 but still 2.1 kcal mol-’ below two non-interacting molecules. The Cl. - -Cl distance is 0.373 nm, and the C-Cl. **Cl angles are 8,-- 137’ and t9,= 100’. The HH constellation CX* * *X-C is commonly encountered in crystals [6]. Most frequently, 8, # 0,: when 8, is close to 180’) 0, gets close to 90’. Our prototypes are too simple perhaps to model numerically the complex structures in the crystallographic literature. Yet they converge persistently on unequal 8 values in the correct ranges. As further examples, we get 151’ and 69’ for HH-MeCl, 130” and 101’ for HH-Me,CHCl. What is it that drives 8, and 0, to differ? This may have to do with the drive towards surface-complementarity [ 211. Traditionally, one has looked for complementarity of shape [ 221, checking whether “the projections of one molecule fit into the hollows of another” [ 231. Nowadays, two other propensities have been identified in the interactions of complex molecules, to wit, electrostatic complementarity and matching of non-polar regions. In many cases, simultaneous fitting of shape and of charges cannot be met. If electrostatic complementarity is imposed (Figs. 1,2), the surfaces of the partners will not lock. If, however, one exploits the translational barrier to prevent electrostatic matching, the computational process converges on the second choice. The optimized channel is then of uniform width (Fig. 4). l
COMPONENTS OF ENERGY-GAIN
In algorithms devolving from Allinger’s MM [ 81, the computed energy comprises six components, due to: (a) bond-length variation; (b) bond-angle variation; (c ) non-bonded interaction (excluding electrostatic) between atoms separated by three bonds; (d) non-bonded interaction (excluding electrostatic) between other atoms; (e) torsion; (f) electrostatic interaction. Denote by E0.i the energy calculated for a single molecule i, and by E the energy calculated for an aggregate. The gain dE is E-C(i)j& (for a pair, i = 1,2 ) , and the well-depth is E= -LIE. For an aggregate, it is convenient to repartition LIE as follows. Skeletal energy (E,, components a, b, c, e). LIE, is strictly intramolecular,
88
always positive and small. It is an indirect measure of the geometrical response of the partners to the imposed proximity. Dispersion energy (Ed, component d). AEd is essentially intermolecular, always negative, and the leading term in dE. It is meant to comprise all interpartner interactions, save those explicitly labeled as “electrostatic”. Electrostatic energy (E,,, component f). LIE,, is essentially intermolecular, quite small for uncharged organic molecules, negative or positive. Table 1 lists e-values for some molecules. There is a way to check their order of magnitude. A rough rule relates the energy of vaporization (L) of a liquid to the well-depth and the number of neighbors (z) closest to a specific molecule: L- 2e/2 [24a]. By using E from computation and L from measurement [ 251, the apparent number of neighbors, z,~~N2 (L/E), is obtained. One may then check whether it falls in a reasonable range. However, the formula cited cannot be used directly. First, L is affected by second-neighbor interactions. To correct for these, E has to be multiplied by a factor > 1 (1.2 N 1.3 [24b,26] ). Second, not all nearest-neighbors can have the same orientation. If Erefers to the deepest minimum, it should be multiplied by a factor < 1. Third, because of thermal vibration at the temperature to which L alludes, the interacting molecules are not at the well-trough, and the multiplier is again < 1. A very rough estimate of this factor comes from the physics of the classical harmonic oscillator [ 271. Here, the mean kinetic energy, i.e., height above the bottom, is half of the total energy (W/2 out of kT). Hence, the factor is - 0.5. An overall multiplier of 3/2 was settled on, viz. zaPP= 3 (L/c). This somewhat arbitrary choice does not affect the check. The values we get for z,,,~ (4 - 10, Table l), do fall in the anticipated range. Also, their evolution is reasonable. Thus, globularity enhances aggregation (cf. pentane and neopentane), and replacement of a light by a heavy halogen deepens the well. TABLE 1 Well depths”
1 Butane 2 Pentane 3 Neopentane 4 Cyclohexane 5 MeCl 6 Me&F-A, 7 Me&Cl-A, 8 Me&Br-A2 9 Me&I-A2
L
E
5.8 6.6 5.6 7.8 5.4
3.4 4.5 3.3 3.9 1.7 3.0 3.5 3.6 4.0
6.9
Z*w 5.1 4.4 5.1 6.0 9.5 5.9
“L, heat of vaporization, kcal mol-‘, from ref. 25; c, calculated well depth, kcal (mol of pairs) -I; Zapp,apparent number of nearest neighbors.
89
To illustrate the relative significance of dE components, here is the breakdown for (Me&Cl) ,-constellations:
A& A& A&,
A2
Al
HT
HH
0.01 -3.28 -0.20
0.09 - 2.40 -0.15
0.02 - 1.57 -0.14
0.00 -2.31 +0.17
Conformational mobility adds a dimension to the potential surface. The simplest case is butane, where liquefaction is known to be attended by coiling [ 14,28,29]. There is evidence that the phenomenon is general to polymethylene chains [ 91. The 11most evident constellations of pairs of butane molecules were optimized, anti (a) and gauche (g), in the combinations au, ag and gg. Data on the calculated energies, assembled in Table 2, may be summarized as follows. (a) All combinations are stabilized by association, but there are energy zones. aaPairs are lowest (relative energy range O-1.2) and then come gu (1.2-2.2) and gg (2.3-3.2). If so, the observed chain-kinking is essentially a mixing effect, due to the larger number of conceivable gu-minima. (b ) It is more likely to find a molecule of g-butane adjacent to u-butane than to another molecule of gbutane. g-Molecules do not tend to cluster. TABLE 2 Calculated energy of butane dimers Energy”
Energy gainb
UUl aa oa3
0.73 1.16 0
- 3.45 -3.02 -4.18
go1 go2 go3 go4 go5
2.15 2.05 1.18 1.61 1.99
-2.67 -2.77 -3.64 -3.21 - 2.83
g&Y1 g@ gg3
2.35 3.15 2.91
-3.11 -2.32 - 2.55
Energies in kcal (mol of pairs ) -I. “Relative to the calculated energy (4.14) of the lowest constellation encountered (~23). bWith respect to the calculated energy of the appropriate non-interacting units.
The components of LIE (energy gain) in the best constellations of each type are as follows:
AE. A& A&.
aa
m3
&?l
+ 0.08 -4.27 +0.01
+0.24 -3.88 -0.00
+0.07 -3.17 -0.01
Here, and in all other cases, the leading contributor is intermolecular attraction. As with single molecules [ 301,some skeletal distortion (reflected in AI?,) mitigates strain and is beneficial to the structure as a whole. Having described results for a prototype (Me&X) with a slight electrostatic bias and for a prototype (C,H.,,) with conformational diversity, we turn to MeCHClCHClMe, where both factors obtain. The anti-gauche equilibrium in meso-2,3dichlorobutane is well characterized. There is about 90% of anti in the vapor, 70% in CC& solution, and 50% in the neat liquid [31,32 1. Calculations by the model of the continuous medium [33] lead to the following interpretation. On condensation, dipolar effects stabilize the gauche form with respect to the anti, and quadrupolar effects destabilize it. Since the former are about twice as ponderant as the latter (1.55 vs. -0.84 kcal mol-’ in the neat liquid), the equilibrium guuchez%nti shifts to the left. Calculations were performed by starting with an optimized butane-butane pair, substituting four hydrogens by chlorines as appropriate, and re-optimizing. First a few gu-pairs with antiparallel C-Cl bonds were examined, derived from the lower energy gu-pairs of butane. The best case corresponded to an energy-gain IZof 4.61 kcal mol-‘. Next, a scan was made of au-constellations, in a search for a case with 6 > 4.61 kcal mol-l: none was encountered. The best au-dimer has ~~4.37 kcal mol-‘. With due reservations, then, the results are consistent with the known increase in g-form on going from vapor to liquids. The surprise lies with the components of LIE. Table 3 presents the breakdown of the calculated energies E and of the energy gains for the two lowest constellations of each type. The parent butane pair, as labeled in Table 2, is given in parentheses. One sees that, when the effects of aggregation are studied by molecular mechanics proper, the calculated electrostatic factor LIE,, is small or negligible. This contrasts with the habitual approach of molecular mechanists to medium effects [ 32,331, where all factors other than electrostatic are assumed negligible. True, in developing a force field, the investigator has a certain latitude of choice as regards the weight of the electrostatic term [ 33,341. The field used here subdues this term by using an intra-body dielectric constant (D) of 4 [ 331. As the numbers in Table 3 show, E,, is significant but LIE,, is not. MM2, with D = 1.5 [ 351, would have dE,,-values some 2-3 times higher than those in the Table, with concomitant small changes in dEd. Even so, AEd
91 TABLE 3 Energy components for the lowest constellations of (MeCHClCHClMe), aa(aa1)
aa(ad)
@a?~)
ga (w4 1
Ed E, Total
10.12 -4.81 - 1.85 3.46
10.58 - 5.56 - 1.87 3.15
11.49 -5.47 - 1.30 4.72
11.40 -5.63 - 1.50 4.27
A& A& A&, Total
0.25 -4.28 -0.03 -4.06
0.71 -5.02 -0.06 -4.37
0.26 - 4.52 +0.11 -4.15
0.17 -4.68 -0.10 -4.61
E
Energies in kcal (mol of pairs)-‘. First 4 lines: as calculated; next 4 lines: relative to the energy of the appropriate non-interacting units.
would lead by far. Moreover, non-bonded interactions have been parametrized such that an electrostatic term is not essential for organic halogen compounds [ 361. Such fields have dE,, z 0 and even E,, E 0. In brief, the two “complementary” tools of the molecular mechanist contradict each other. While the numerical demonstration stresses A&, it does not affect the conceptual notion of “electrostatic complementarity”. Geometries of the best aggregates are such that electronegative centers approach alkyl residues. Headto-head orientations are stabilized if C-X*. lX angles can be obtuse. The reason probably lies with propensities reflected in the Hill equation [ 341. Efficient attraction requires both a deep well and the proximity of interacting groups. Numerical details are such that C-X. - *H-C is preferable to C-X. **XC and to C-H. - *H-C. SIZE AND SHAPE
When two molecules approach each other to form a stable aggregate, the geometrical parameters of each (bond lengths, angles) vary somewhat [5]. Two types of geometrical descriptor are therefore of relevance: those that describe the aggregate as a whole, and those that describe changes in each separate partner. For the first type, descriptors originally developed to characterize the van der Waals’ body of single molecules [ 41 are appropriate. These are the volume of the body ( VW), the surface area (S,) , the radius r, of a sphere of volume equal to VW [ 371, the globularity 4me2/S, [ 2,4]. In Fig. 5, a circle of radius r, is superposed on the body of meso-anti-MeCHClCHClMe. Since molecules are better emulated by ellipsoids than by spheres, it should be advantageous to refer to an “equivalent ellipsoid”. This is an ellipsoid of revolution, of volume VW,centered at the molecular center, and exactly circumscribing the atom farthermost from the center. The circumscription constraint
92
Fig. 5. anti-meso-2,3-Dichlorobutane. Figs. 5-7 show the equivalent sphere, the generating ellipse (internal), an ellipse removed from it by 0.3 nm (external), and the non-bonded energy (thickened curve) of an H-probe revolving along the latter.
fixes the length and direction of the major axis of the generating ellipse, and the volume constraint fixes the length of the minor axis b= (3Vw/4na)‘/2 The inner ellipse in Fig. 5 is the generating ellipse of the “equivalent ellipsoid”. Its eccentricity provides an immediate apprehension of molecular shape. When expressed numerically e= [l-
(b/a)2]‘/2,
it becomes a descriptor of a type required in the electrostatic theory of liquids [ 381. Also, the surface area of the ellipsoid is easily calculated [ 391 &,,=2n(b/a)(ab+
(a3/(a2-b”)““)
arcsin( (a2-b2)l12/a)).
Since the surface of any molecule is somewhat convoluted, S,,, is habitually smaller than S,. The protrusions of the molecular body through the ellipsoid surface provide for visual apprehension of surface convolution. When expressed numerically, the ratio
becomes a descriptor of surface roughness, required in the theory of surface accessibility [ 401. Figures 6 and 7 show, respectively, the best aa-dimer and the best ga-dimer of MeCHClCHClMe. On going from the monomer to the dimers, eccentricity increases (e changes from 0.777 to, respectively, 0.819 and 0.860), as does the surface convolution (SJS,,, changes from 1.215 to, respectively, 1.510 and 1.475). A counter example, where eccentricity decreases on going from the
UNIT:
0.1
NM OR
2.5
KJ/MOL
Fig. 6. Best aa-dimer of MeCHClCHClMe. Fig. 7. Best gu-dimer of MeCHClCHClMe.
L-_--X UNIT: 0.1
NM
Fig. 8. Parallel alignment of two all-anti pentanes. Fig. 9. anti-meso-2,3-Dichlorobutane, revolving along the external ellipse.
showing the electrostatic energy of a positive unit of charge,
monomer to the dimer, is provided by pentane. In the all-anti conformation of pentane, e = 0.880 and S,/S,ll = 1.153. For the dimer shown in Fig. 8, e=0.847 and SW/Sell= 1.489. In Fig. 8, note the uniformity of the intermolecular channel, i.e., the way “bumps pack in hollows”. Since calculations stress non-bonded intermolecular interactions (AEd, Table 3)) it was of interest to compute the non-bonded energy of an atomic probe revolving beyond the generating ellipse. An hydrogen atom was chosen as probe, and Hill-equation parameters were as in MM2 [41]. The chosen trajectory was an ellipse (external ellipse in Figs. 5-7), defined by producing the semiaxes of the generating ellipse by a distance 1.At Z= 0.2 nm, the energy is highly repulsive. At 0.25 nm, it is mildly attractive. In Figs. 5-7, the external ellipse
94
is drawn at 1=0.3 nm, less than the order of intermolecular internuclear distances in the dimers (closest approaches N 0.4 nm). The computed energy on the trajectory is traced as a thick curve [ 421, where the trajectory serves as the coordinate zero. Thus, points within the external ellipse (all points in Figs 57, where 1= 0.3 nm) correspond to attraction. Unlike tracings of the molecular electrostatic potential (MEP), the curve obtained for the non-bonded energy is unstructured: Fig. 5 cannot predict the constellations in Figs. 6 and 7. One cannot claim that, in this context, shape is of minor significance [ 43 1. Rather, unlike in computing MEPs, a probe cannot be representative of an actual environment in a neighboring molecule. On replacing the H-probe by, say, a Cl-atom in a second molecule, two parameters change in the Hill equation (r* and E [8,34] ) and the computed energy is not linear in either. Finally, one is curious to know why electrostatics play such a minor role. To this end we inspected the electrostatic energy of a revolving unit of positive charge. In Fig. 9, this is done for the molecule of Fig. 5, again at a distance of 0.3 nm from the generating ellipse. The charge distribution has been estimated by the modified Del Re method [44]. The ensuing curve does have structure, and minima (arrowheads attracting the probe) are well defined. Maxima (humps repelling the probe) are however weak and shallow. Unlike the nonbonded energy, the electrostatic energy is proportional to the energy of the probe. For an actual atom in a neighboring molecule, charged as it actually is (Qcl,- 0.168; qn, 0.4 - 0.6 electron units), all fine structure would level off. INTRAMOLECULAR
CHANGES
Now to descriptors of the second type. The question to be answered is as follows: when brought to a given environment, does a given molecule contract or expand? Again we consider the van der Waals’ body. Atoms closer to each other than the sum of assigned radii form an overlap lens. Its volume is a measure of their proximity. It is meaningful to distinguish two types of lens: (a) formed by bonded atoms and reflecting bond-lengths; (b ) formed by non-bonded atoms and reflecting bond angles or dihedral angles. Usage of van der Waals’ radii is advantageous, since atoms removed from each other farther than A and D in a bonded sequence A-B-C-D do not overlap. Let On be the sum of bonded overlaps and ONBthe sum of non-bonded overlaps. In analogy with the definition of energy gains, define now “overlap gains”: don=o,-c(i)oO,,~ dONn=ONn-Z(i)O”Nn,i (for bimolecular association, i = 1,2). Positive values express an overall con-
95
TABLE 4 Boeyens’ analysis for compounds of stereochemical interest 0,
W-L
se
0 NB
E
ccl CaH&le, eq ax
39.600 132.499 132.483
9.219 9.142 53.696 53.654
- 640.5 - 639.7 - 1898.8 - 1896.9
Decalin, trans
187.532 187.503
83.691 83.506
- 2615.3 -2612.3
43.208 43.187
10.892 10.895
- 605.6 - 605.0
39.623
CiS
ClCH&H$l,
a iz
Overlaps of van der Waals’ spheres in nm3 x 10e3 (i.e., A 3); energies of Boeyens’ overlaps in kcal mol-‘.
traction of the associated components; negative values express an overall expansion. Atomic overlaps were studied by Boeyens, who also used them to evaluate bond energies [ 45,461. Clearly, the usefulness of the tool goes beyond the few computational examples here at hand. To put numbers into perspective, some results of stereochemical interest are given in Table 4. All refer to computed geometries (alkanes by MM2, chlorides as in ref. 3). The overlaps are by van der Waals radii [lo], and the total energies E by Boeyens radii [ 451. On passing from a less to a more congested stereoisomer, dOs and LION are negative, reflecting expansion. This is because bonds generally stretch and valence angles open up. Actual numbers are small. In the two decalins, for example, d (0s + ONB) is only 0.214 x 10m3nm3 molecule-‘. This is because interatomic distances change only slightly when compared with atomic radii, and the overlap volume depends only weakly on the interatomic distance. In fact, the overlap formula [45] leads to [d(overlap)]/[d(distance)]
=n (; d+r)(+)
where d is the distance, and r is the radius assigned to the overlapping atoms. Now referring to the dimers, in all cases studied (butane, pentane, neopentane, cyclohexane, MeCHClMe, Me3CX, MeCHClCHClMe ) , computed geometries correspond to contraction. It is as if intramolecular squeezing enhances the chance of intermolecular attraction. The numbers, again, are small. For the best au-constellation of (MeCHClCHClMe), (Fig. 6), d(On+ONB)=0.207X10-3nm3molecule-’of the same order as the decalin contraction. The corresponding number for the best gu-constellation, presumably the global minimum (Fig. 7)) is even lower (0.109 x 10S3 nm3 molecule- ’ ) .
dE, also indicates a higher distortion in au (Table 3)) but cannot pinpoint the sense of distortion. CONCLUSION
It has recently been stated that “weak non-covalent bonds involving neutral organic molecules are of vital importance, e.g., in molecular biology, drug design, etc.” [ 471. The statement underlies a curious fashion. Computationalists tend to study molecular association either when the partners are quite complex (e.g., &trypsin with pancreatic trypsin inhibitor [ 21b] ), or when the partners are quite simple (e.g., NH:, with HF, HCI with Nz, NH3 with LiF [48] ). The few examples given above serve to suggest that examination of molecules in between is not unrewarding. In particular, it brings to the fore questions of principle that run the risk of being ignored in dealing with the extremes of the complexity scale. Very small dimers can be treated effectively by quantum chemistry [ 481. Molecular mechanics, though created specifically to meet single molecules, can be applied to more complex dimers. No attempt was made here to optimize the force field, and improvements can be envisaged. As it stands, it did characterize head-to-head orientations of organic halogen compounds, the intramolecular response to intermolecular interactions, and the kinking of carbon-chains in condensed phases. As for descriptors of size and shape, it is true that their content of information becomes useful mainly in the context of complex molecules [491. But it is through application to molecules of lower complexity that they can be tested to advantage [ 41. Here we considered descriptors that reflect the overall change in geometry, eccentricity and surface roughness, on going from a single molecule to an associated pair.
REFERENCES 1 2 3 4 5 6 7 8 9 10 11
A.Y. Meyer, J. Chem. Sot. Perkin Trans., 2 (1986) 1567. A.Y. Meyer, J. Comput. Chem., 9 (1988) 18. A.Y. Meyer, J. Mol. Struct., 94 (1983) 95, and previous papers in the series. A.Y. Meyer, Chem. Sot. Rev., 15 (1986) 449. I. Hargittai and M. Hargittai, in Molecular Structures and Energetics, Vol. 2, VCH, Weinheim, 1987. N. Ramasubbu, R. Parthasarathy and P. Murray-Rust, J. Am. Chem. Sot., 108 (1986) 4308. N.L. Allinger, M.A. Miller, F.A. Van Catledge and J.A. Hirsch, J. Am. Chem. Sot., 89 (1967) 4345. N.L. Allinger,M.T. Tribble, M.A. Mil1erandD.H. Wertz, J.Am. Chem. Soc.,93 (1971) 1637. J. Gajdoli and T. Bleha, Collect. Czech. Chem. Commun., 48 (1983) 71. A. Gavezzotti, J. Am. Chem. Sot., 105 (1983) 5220. J.O. Hirschfelder, J. Chem. Phys., 8 (1940) 431.
97 12 13 14 15
16
17 18 19 20 21
22 23 24 25 26 27
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
J.L. Letelier, Y.-N. Chiu and G. Zon, J. Mol. Struct., 63 (1980) 273. U. Burkert and N.L. Aihnger, J. Comput. Chem., 3 (1982) 40. K.B. Wiberg and R.H. Boyd, J. Am. Chem. Sot., 94 (1972) 8426. Since the terms “conformation” and “configuration” have specific meanings in the context ofsingle molecules, we adopt Prelog’s “con&elk&ion” when referring to the mutual orientation of two molecules. See V. Prelog, J. Chem. Sot., (1950) 420. This exotic energy unit was chosen because enthalpies of vaporization, relevant to the discussion and for comparison, are traditionally tabulated in kcaI mol-‘. Excepting its first occurrence and in tables, the unit is abbreviated as kcai mol-‘. E. Hirota, Y. Endo, S. Saito and J.L. Duncan, J. Mol. Spectrosc., 89 (1981) 285. A. Bondi, J. Phys. Chem., 68 (1964) 441. V. Cody and P. Murray-Rust, J. Mol. Struct., 112 (1984) 189. 0. Bastiansen, K. Kveseth and H. MsllendaI, Topics CUR. Chem., 81 (1979) 99. Complementarity has been studied extensively by G. Naray-Sxab& See, for example, (a) G. Nbray-SxabQ J. Am. Chem. Sot., 106 (1984) 4584; (b) G. Naray-Szab6 and P. Nagy, Enzyme, 36 (1986) 44; (c) P. Nagy, J.G. Angyan and G. Naray-Sxab6, J. Mol. Struct., 149 (1987) 169. M. Santavjr and J. Kypr, J. Mol. Graphics, 2 (1984) 47. A.I. Kitaigorodski, Chem. Sot. Rev., 7 (1978) 133. D. Tabor, Gases, Liquids and Solids, 2nd edn., Cambridge University Press, Cambridge, 1979. (a) p. 153; (b) p. 177. R.C. Weast (editor), Handbook of Chemistry and Physics, 52nd edn., The Chemical Rubber Co., Cleveland, OH, 1971-1972,~. D151. T.L. Hill, Lectures on Matter and Equilibrium, Benjamin, New York, 1966, p. 54. F. Reif, Fundamentals of Statistical and Thermal Physics, McGraw-Hill, New York, 1965, p. 251. I am indebted to Professor S. Baer, of the Hebrew University, for bringing this analysis to my attention. L.R. Pratt, C.S. Hsu and D. Chandler, J. Chem. Phys., 68 (1978) 4202. P.E. Schoen, R.G. Priest, J.P. Sheridan and J.P. Schnur, J. Chem. Phys., 71 (1979) 317. F.H. Westheimer, in MS. Newman (Ed.), Steric Effects in Organic Chemistry, Wiley, New York, 1956, pp. 523-555, in particular p. 550. A.A. Bothner-By and C. Naar-Colin, J. Am. Chem. Sot., 84 (1962) 743. R.J. Abraham and E. Bretschneider, in W.J. Orville-Thomas (Ed.), Internal Rotation in Molecules, Wiley, London, 1974. A.Y. Meyer, J. Comput. Chem., 2 (1981) 384. A.Y. Meyer and F.R.F. Forrest, J. Comput. Chem., 6 (1985) 1. A.Y. Meyer, N.L. Allinger and Y. Yuh, Israel J. Chem., 20 (1980) 57. T. Rydiand and R. Stralevik, J. Mol. Struct., 105 (1983) 157; R. Stelevik, J. Mol. Struct., 109 (1984) 397. J.T. Edward and P.G. Farrell, Can. J. Chem., 53 (1975) 2965. A. Raudino, F. Zuccarello and G. Buemi, J. Chem. Sot., Faraday Trans. 2,77 (1981) 1331. See, for example, W.L. Duren, Calculus and Analytic Geometry. Xerox College publishing, Lexington, MA, 1972, p. 700. J.L. Finney, J. Mol. Biol., 119 (1978) 415. N.L. Allinger, Operating Instructions for MM2 and MMP2 Programs, updated as of August 1985, University of Georgia, 1985. Curve-thickening was achieved by the routines of P. Senn, Comp. Chem., 10 (1986) 219. J.S. Rowlinson, Discuss. Faraday Sot., 49 (1970) 30. A.Y. Meyer, J. Chem. Sot., Perkin Trans., 2 (1982) 1199. J.C.A. Boeyens, J. South Afr. Chem. Inst., 26 (1973) 94. J.C.A. Boeyens, South Afr. J. Chem., 33 (1980) 14, and previous publications in the series.
47 48
49
E. Weber, in the Preface to Top. Curr. Chem., Vol. 140, Springer-Verlag, Berlin, 1987. For very recent examples of ab initio work on dimers, see: G.J.B. Hurst, P.W. Fowler, A.J. Stone and A.D. Buckingham, Int. J. Quantum Chem., 29 (1986) 1223; B. Silvi and M. Allawena, J. Mol. Struct., 135 (1986) 225; Z. Latajka, K. Morokuma, H. Ratajczak and W.J. Orville-Thomas, J. Mol. Struct., 135 (1986) 429. D.E. Walters and A.J. Hopfinger, J. Mol. Struct., 134 (1986) 317.
MOLECULAR
MECHANICS AND MOLECULAR
83
SHAPE
Part VI. The response of simple molecules to bimolecular association****
A.Y. MEYER Department of Organic Chemistry, Hebrew University, Jerusalem (Israel) (Received 20 July 1987)
ABSTRACT A molecular-mechanical force field, originally developed to calculate the structure and conformation of single molecules, has been applied to pairs of identical molecules lying at close proximity. The prototypes examined are very simple - short alkanes and short haloalkanes - but do include cases with conformational diversity and possible electrostatic effects. The optimized mutual orientations include the parallel alignment of linear alkanes, head-to-tail orientations of haloalkanes, as well as the head-to-head association, known from the crystallographic literature. The effect of aggregation on the conformation of the separate partners has been studied, and the intramolecular response characterized: in all cases examined, the partners contract. The computed energy gains are of the expected order of magnitude, the main contribution coming always from what molecular-mechanists label “non-bonded attraction”. Graphical representations of optimized orientations show the “protrusion-in-cleft” matching of molecular surfaces.
INTRODUCTION
The author has reported on the molecular mechanics [3] and on the size and shape [1,2,4] of organic halogen compounds. Two trends in the current literature make it timely to go beyond single molecules, and apply the existing computational techniques to bimolecular aggregates. One is the rising interest in intramolecular geometrical changes that accompany molecular association [5]. The other is the accumulation of data regarding angular preferences of intermolecular forces [ 61. The aim here is not to re-optimize force field parameters and reproduce documented data. On the contrary, it is to examine to what extent a single-molecule force field can, together with appropriate descriptors of size and shape, contend with association and add to the conceptual background. *Dedicated to Professor Bernard Pullman. **For Parts IV and V, see refs. 1 and 2.
0166-1280/88/$03.50
0 1988 Elsevier Science Publishers B.V.
Thus, the questions we ask are primarily of principle, in the spirit of former molecular-mechanical examinations of the hexane crystal [ 7,8] and of hexaneassociates [ 91. To get answers, however, extensive numerical checking is unavoidable. As a start, pairs of very simple molecules were examined, alkanes and haloalkanes. The force field has been described [3]. Use of a different computational technique, molecular-mechanical or another, would undoubtedly affect numerical detail. It should not affect the lines of approach, nor, with evident reservations, the nature of results. The main difference of our force field from others is that all bonds (including C-H) are attributed bond dipoles, except C-C bonds when neither carbon atom carries a heteroatomic substituent. The descriptor-evaluating program has also been described [ 1,2,4]. The literature has accorded little attention to dimers of the types to be discussed. We therefore find it helpful to provide computer-drawn cuts through the van der Waals’ bodies of the dimers. To obtain these, each atom is represented by a sphere of the appropriate van der Waals’ radius [lo], each molecule thus becoming a system of interlocking spheres. To help in viewing, nuclei within 0.04 nm of the plane of cut are marked by circles, of radius equal to half the atomic radius. Nuclei farther away are not shown. Sparse hatching indicates an heteroatom. To convey a feeling for aggregate size, all figures include the “equivalent ellipse”. This is the generating ellipse of an ellipsoid of revolution that (a) exactly circumscribes the atom farthest away from the center, and (b) whose volume is equal to the volume of the van der Waals body. In most figures, the center is center of mass, and the plane shown is the first inertial plane (i.e., defined by 1, and 1s [ 111, here denoted xy). In a few instances we preferred the center of volume and the first plane by volume. These are obtained exactly as in the inertial case, except that atoms are weighted by cubed van der Waals radii. EXTREMA ON THE POTENTIAL SURFACE
As is well known, in “optimizing a chemical structure”, the investigator starts by guessing a “trial geometry”. Next, the computer program is allowed to modify the geometry by steps, such that energies computed for sequential structures go steadily down. The process stops at a “final geometry”, whereat further lowering is not feasible. Since minima are separated by ridges and the process works downhill, a minimum risks being overlooked if the trial geometry is not within the proper valley. In dealing with internal rotation in single molecules, constraints of chemical bonding, previous knowledge, and symmetry considerations [ 121 serve as guides. When in doubt, one may “drive” along the saddles [ 13,141. By contrast, little is known about preferred relative orientations of two (or more) molecules, and detailed scans of potential surfaces - allowing relaxation in all but a few degrees of freedom - are not worth the while. Consider the constellations [ 151 of Me&X-dimers (X = halogen). Me&X
is conformationally uniform and highly symmetric, and has electrostatically defined head (X) and tail (Me&). One expects a minimum when the C-X moieties of the two molecules are antiparallel to each other. Computation reveals splitting in two separate minima. One is “singly flanking”, A,: the Clatom in one sole partner is flanked by two methyl groups of the other partner (Fig. 1). The other is “doubly flanking”, AZ: the Cl-atom in both partners is thus flanked (Fig. 2). AZ is a deeper minimum and, for all we know, the minimum minimorum. Computed energy differences, E (A, ) -E (A,), are 0.8 kcal (molofpairs)-’ [16] forX=F,1.0forX=C1,Br,and0.3forX=I.A1andA, are separated by a barrier to external rotation. When a residue of local symmetry C,, rotates internally (e.g., internal rotation in Me-Me), the periodicity is 3, there are 3 minima in a cycle and 6 extrema in all. In the Fourier series V(m) = ( V,/2) (1 -cos 3~) + ( V,/ 2) (l-cos 6w)+ .**, the term in V, has the correct symmetry and is leading [ 171. By contrast, the external rotation A,*A,* ** has 6 minima in a cycle and 12 extrema in all. The term in V, does not reflect correctly the sequence of events, and calculations suggest that it is not leading. More specifically, assume that maxima occur halfway between A, and A2 (rotated through 46 with respect to either), and denote by dA and dB, respectively, the energy of Al and of the barrier above AZ. In the Fourier series, V, =dA and V, =dB - (dA/2). To estimate dB in Me&I, one of the molecules in constellation A, was rotated through 46 about C-I. The energy rose to Y
Y
L
-----
UNIT: 0.1
NM
XL_--__X UNIT: 0.1
NM
Fig. 1. “Singly-flanking” constellation (A,) of (Me&I),. In all figures the molecule is cut across the first inertial plane, except in Figs. 5 and 9 where this is the first plane by volume. All figures include the generating ellipse of the equivalent ellipsoid. Fig. 2. “Doubly-flanking”constellation
(A,) of (Me&I),.
86
2.57 kcal mol-l. Holding in place the central carbon atom and one other carbon in both units, the assembly was allowed to relax. dB fell to 0.98 kcal mol-l. If the higher estimate is the more realistic [ 81, V, - 0.3 and V, - 2. At any rate, V, is very much greater than V,. The geometry at the barrier is shown in Fig. 3. In all three constellations (Figs. l-3), the van der Waals’ body of a Me&I-unit has virtually the same size: its volume is N 0.109 nm3 molecule-‘, its surface area is N 1.35 nm2 molecule-i. What the constellation does affect are the dimensions of the intermolecular channel. This is the “bimolecular variant” of a conclusion drawn for single molecules: configurational isomers differ in density mainly because of differences in packing, less so because of differences in size [ 41. At the barrier (Fig. 3), the partners touch or interpenetrate, and the channel is blocked. It should be stressed that the energy at the barrier is still well below that of two non-interacting molecules. At least in this case, external rotation within the aggregate is cheaper to activate than dissociation followed by reorientation. Generalizing now, it may happen that normal contacts along extensive ranges of the interacting surfaces require a close contact at some point. Then, it may again happen that stabilization due to the former outweighs the destabilization caused by the latter. Now, the van der Waals’ volume of a molecule was defined in 1964 as the volume that is impenetrable to other molecules with thermal energies at ordinary temperatures [ 181. Later work does not substantiate this view. Recent scans of the Cambridge Crystallographic Data Base revealed 729 43 =
l______x UNIT:
0.1
NM
L-_-_-X UNIT:
0.1
Fig. 3. Barrier separating constellations A, and A2 of (Me,CI),. Fig. 4. “Head-to-head” constellation
(HH) of (Me&Cl),.
NM
instances of interpenetrating C-X. * (0,N) contacts [6] and 86 instances of interpenetrating C-1. **(O,N,S) contacts [ 191. An intramolecular instance is provided by 2,2 ’ -dihalobiphenyls, where halogen atoms approach closer than the sum of theirvan der Waals’s radii [ 201. Our simple prototype shows that intdr!enetration does not imply an energy-rise. It may well enable a considerable energy-gain. Two other minima were located on the potential surface of (Me&X)2, this time separated from A, (Fig. 2) by barriers to translation. In one, the partners are head-to-tail (HT, Me&X* - *Me&X); in the other, they are head-to-head (HH, Me&X* * *XCMe,). The more interesting constellation, HH, is shown in Fig. 4. By computation, it lies 1.3 kcal mol-’ above A2 but still 2.1 kcal mol-’ below two non-interacting molecules. The Cl. - -Cl distance is 0.373 nm, and the C-Cl. **Cl angles are 8,-- 137’ and t9,= 100’. The HH constellation CX* * *X-C is commonly encountered in crystals [6]. Most frequently, 8, # 0,: when 8, is close to 180’) 0, gets close to 90’. Our prototypes are too simple perhaps to model numerically the complex structures in the crystallographic literature. Yet they converge persistently on unequal 8 values in the correct ranges. As further examples, we get 151’ and 69’ for HH-MeCl, 130” and 101’ for HH-Me,CHCl. What is it that drives 8, and 0, to differ? This may have to do with the drive towards surface-complementarity [ 211. Traditionally, one has looked for complementarity of shape [ 221, checking whether “the projections of one molecule fit into the hollows of another” [ 231. Nowadays, two other propensities have been identified in the interactions of complex molecules, to wit, electrostatic complementarity and matching of non-polar regions. In many cases, simultaneous fitting of shape and of charges cannot be met. If electrostatic complementarity is imposed (Figs. 1,2), the surfaces of the partners will not lock. If, however, one exploits the translational barrier to prevent electrostatic matching, the computational process converges on the second choice. The optimized channel is then of uniform width (Fig. 4). l
COMPONENTS OF ENERGY-GAIN
In algorithms devolving from Allinger’s MM [ 81, the computed energy comprises six components, due to: (a) bond-length variation; (b) bond-angle variation; (c ) non-bonded interaction (excluding electrostatic) between atoms separated by three bonds; (d) non-bonded interaction (excluding electrostatic) between other atoms; (e) torsion; (f) electrostatic interaction. Denote by E0.i the energy calculated for a single molecule i, and by E the energy calculated for an aggregate. The gain dE is E-C(i)j& (for a pair, i = 1,2 ) , and the well-depth is E= -LIE. For an aggregate, it is convenient to repartition LIE as follows. Skeletal energy (E,, components a, b, c, e). LIE, is strictly intramolecular,
88
always positive and small. It is an indirect measure of the geometrical response of the partners to the imposed proximity. Dispersion energy (Ed, component d). AEd is essentially intermolecular, always negative, and the leading term in dE. It is meant to comprise all interpartner interactions, save those explicitly labeled as “electrostatic”. Electrostatic energy (E,,, component f). LIE,, is essentially intermolecular, quite small for uncharged organic molecules, negative or positive. Table 1 lists e-values for some molecules. There is a way to check their order of magnitude. A rough rule relates the energy of vaporization (L) of a liquid to the well-depth and the number of neighbors (z) closest to a specific molecule: L- 2e/2 [24a]. By using E from computation and L from measurement [ 251, the apparent number of neighbors, z,~~N2 (L/E), is obtained. One may then check whether it falls in a reasonable range. However, the formula cited cannot be used directly. First, L is affected by second-neighbor interactions. To correct for these, E has to be multiplied by a factor > 1 (1.2 N 1.3 [24b,26] ). Second, not all nearest-neighbors can have the same orientation. If Erefers to the deepest minimum, it should be multiplied by a factor < 1. Third, because of thermal vibration at the temperature to which L alludes, the interacting molecules are not at the well-trough, and the multiplier is again < 1. A very rough estimate of this factor comes from the physics of the classical harmonic oscillator [ 271. Here, the mean kinetic energy, i.e., height above the bottom, is half of the total energy (W/2 out of kT). Hence, the factor is - 0.5. An overall multiplier of 3/2 was settled on, viz. zaPP= 3 (L/c). This somewhat arbitrary choice does not affect the check. The values we get for z,,,~ (4 - 10, Table l), do fall in the anticipated range. Also, their evolution is reasonable. Thus, globularity enhances aggregation (cf. pentane and neopentane), and replacement of a light by a heavy halogen deepens the well. TABLE 1 Well depths”
1 Butane 2 Pentane 3 Neopentane 4 Cyclohexane 5 MeCl 6 Me&F-A, 7 Me&Cl-A, 8 Me&Br-A2 9 Me&I-A2
L
E
5.8 6.6 5.6 7.8 5.4
3.4 4.5 3.3 3.9 1.7 3.0 3.5 3.6 4.0
6.9
Z*w 5.1 4.4 5.1 6.0 9.5 5.9
“L, heat of vaporization, kcal mol-‘, from ref. 25; c, calculated well depth, kcal (mol of pairs) -I; Zapp,apparent number of nearest neighbors.
89
To illustrate the relative significance of dE components, here is the breakdown for (Me&Cl) ,-constellations:
A& A& A&,
A2
Al
HT
HH
0.01 -3.28 -0.20
0.09 - 2.40 -0.15
0.02 - 1.57 -0.14
0.00 -2.31 +0.17
Conformational mobility adds a dimension to the potential surface. The simplest case is butane, where liquefaction is known to be attended by coiling [ 14,28,29]. There is evidence that the phenomenon is general to polymethylene chains [ 91. The 11most evident constellations of pairs of butane molecules were optimized, anti (a) and gauche (g), in the combinations au, ag and gg. Data on the calculated energies, assembled in Table 2, may be summarized as follows. (a) All combinations are stabilized by association, but there are energy zones. aaPairs are lowest (relative energy range O-1.2) and then come gu (1.2-2.2) and gg (2.3-3.2). If so, the observed chain-kinking is essentially a mixing effect, due to the larger number of conceivable gu-minima. (b ) It is more likely to find a molecule of g-butane adjacent to u-butane than to another molecule of gbutane. g-Molecules do not tend to cluster. TABLE 2 Calculated energy of butane dimers Energy”
Energy gainb
UUl aa oa3
0.73 1.16 0
- 3.45 -3.02 -4.18
go1 go2 go3 go4 go5
2.15 2.05 1.18 1.61 1.99
-2.67 -2.77 -3.64 -3.21 - 2.83
g&Y1 g@ gg3
2.35 3.15 2.91
-3.11 -2.32 - 2.55
Energies in kcal (mol of pairs ) -I. “Relative to the calculated energy (4.14) of the lowest constellation encountered (~23). bWith respect to the calculated energy of the appropriate non-interacting units.
The components of LIE (energy gain) in the best constellations of each type are as follows:
AE. A& A&.
aa
m3
&?l
+ 0.08 -4.27 +0.01
+0.24 -3.88 -0.00
+0.07 -3.17 -0.01
Here, and in all other cases, the leading contributor is intermolecular attraction. As with single molecules [ 301,some skeletal distortion (reflected in AI?,) mitigates strain and is beneficial to the structure as a whole. Having described results for a prototype (Me&X) with a slight electrostatic bias and for a prototype (C,H.,,) with conformational diversity, we turn to MeCHClCHClMe, where both factors obtain. The anti-gauche equilibrium in meso-2,3dichlorobutane is well characterized. There is about 90% of anti in the vapor, 70% in CC& solution, and 50% in the neat liquid [31,32 1. Calculations by the model of the continuous medium [33] lead to the following interpretation. On condensation, dipolar effects stabilize the gauche form with respect to the anti, and quadrupolar effects destabilize it. Since the former are about twice as ponderant as the latter (1.55 vs. -0.84 kcal mol-’ in the neat liquid), the equilibrium guuchez%nti shifts to the left. Calculations were performed by starting with an optimized butane-butane pair, substituting four hydrogens by chlorines as appropriate, and re-optimizing. First a few gu-pairs with antiparallel C-Cl bonds were examined, derived from the lower energy gu-pairs of butane. The best case corresponded to an energy-gain IZof 4.61 kcal mol-‘. Next, a scan was made of au-constellations, in a search for a case with 6 > 4.61 kcal mol-l: none was encountered. The best au-dimer has ~~4.37 kcal mol-‘. With due reservations, then, the results are consistent with the known increase in g-form on going from vapor to liquids. The surprise lies with the components of LIE. Table 3 presents the breakdown of the calculated energies E and of the energy gains for the two lowest constellations of each type. The parent butane pair, as labeled in Table 2, is given in parentheses. One sees that, when the effects of aggregation are studied by molecular mechanics proper, the calculated electrostatic factor LIE,, is small or negligible. This contrasts with the habitual approach of molecular mechanists to medium effects [ 32,331, where all factors other than electrostatic are assumed negligible. True, in developing a force field, the investigator has a certain latitude of choice as regards the weight of the electrostatic term [ 33,341. The field used here subdues this term by using an intra-body dielectric constant (D) of 4 [ 331. As the numbers in Table 3 show, E,, is significant but LIE,, is not. MM2, with D = 1.5 [ 351, would have dE,,-values some 2-3 times higher than those in the Table, with concomitant small changes in dEd. Even so, AEd
91 TABLE 3 Energy components for the lowest constellations of (MeCHClCHClMe), aa(aa1)
aa(ad)
@a?~)
ga (w4 1
Ed E, Total
10.12 -4.81 - 1.85 3.46
10.58 - 5.56 - 1.87 3.15
11.49 -5.47 - 1.30 4.72
11.40 -5.63 - 1.50 4.27
A& A& A&, Total
0.25 -4.28 -0.03 -4.06
0.71 -5.02 -0.06 -4.37
0.26 - 4.52 +0.11 -4.15
0.17 -4.68 -0.10 -4.61
E
Energies in kcal (mol of pairs)-‘. First 4 lines: as calculated; next 4 lines: relative to the energy of the appropriate non-interacting units.
would lead by far. Moreover, non-bonded interactions have been parametrized such that an electrostatic term is not essential for organic halogen compounds [ 361. Such fields have dE,, z 0 and even E,, E 0. In brief, the two “complementary” tools of the molecular mechanist contradict each other. While the numerical demonstration stresses A&, it does not affect the conceptual notion of “electrostatic complementarity”. Geometries of the best aggregates are such that electronegative centers approach alkyl residues. Headto-head orientations are stabilized if C-X*. lX angles can be obtuse. The reason probably lies with propensities reflected in the Hill equation [ 341. Efficient attraction requires both a deep well and the proximity of interacting groups. Numerical details are such that C-X. - *H-C is preferable to C-X. **XC and to C-H. - *H-C. SIZE AND SHAPE
When two molecules approach each other to form a stable aggregate, the geometrical parameters of each (bond lengths, angles) vary somewhat [5]. Two types of geometrical descriptor are therefore of relevance: those that describe the aggregate as a whole, and those that describe changes in each separate partner. For the first type, descriptors originally developed to characterize the van der Waals’ body of single molecules [ 41 are appropriate. These are the volume of the body ( VW), the surface area (S,) , the radius r, of a sphere of volume equal to VW [ 371, the globularity 4me2/S, [ 2,4]. In Fig. 5, a circle of radius r, is superposed on the body of meso-anti-MeCHClCHClMe. Since molecules are better emulated by ellipsoids than by spheres, it should be advantageous to refer to an “equivalent ellipsoid”. This is an ellipsoid of revolution, of volume VW,centered at the molecular center, and exactly circumscribing the atom farthermost from the center. The circumscription constraint
92
Fig. 5. anti-meso-2,3-Dichlorobutane. Figs. 5-7 show the equivalent sphere, the generating ellipse (internal), an ellipse removed from it by 0.3 nm (external), and the non-bonded energy (thickened curve) of an H-probe revolving along the latter.
fixes the length and direction of the major axis of the generating ellipse, and the volume constraint fixes the length of the minor axis b= (3Vw/4na)‘/2 The inner ellipse in Fig. 5 is the generating ellipse of the “equivalent ellipsoid”. Its eccentricity provides an immediate apprehension of molecular shape. When expressed numerically e= [l-
(b/a)2]‘/2,
it becomes a descriptor of a type required in the electrostatic theory of liquids [ 381. Also, the surface area of the ellipsoid is easily calculated [ 391 &,,=2n(b/a)(ab+
(a3/(a2-b”)““)
arcsin( (a2-b2)l12/a)).
Since the surface of any molecule is somewhat convoluted, S,,, is habitually smaller than S,. The protrusions of the molecular body through the ellipsoid surface provide for visual apprehension of surface convolution. When expressed numerically, the ratio
becomes a descriptor of surface roughness, required in the theory of surface accessibility [ 401. Figures 6 and 7 show, respectively, the best aa-dimer and the best ga-dimer of MeCHClCHClMe. On going from the monomer to the dimers, eccentricity increases (e changes from 0.777 to, respectively, 0.819 and 0.860), as does the surface convolution (SJS,,, changes from 1.215 to, respectively, 1.510 and 1.475). A counter example, where eccentricity decreases on going from the
UNIT:
0.1
NM OR
2.5
KJ/MOL
Fig. 6. Best aa-dimer of MeCHClCHClMe. Fig. 7. Best gu-dimer of MeCHClCHClMe.
L-_--X UNIT: 0.1
NM
Fig. 8. Parallel alignment of two all-anti pentanes. Fig. 9. anti-meso-2,3-Dichlorobutane, revolving along the external ellipse.
showing the electrostatic energy of a positive unit of charge,
monomer to the dimer, is provided by pentane. In the all-anti conformation of pentane, e = 0.880 and S,/S,ll = 1.153. For the dimer shown in Fig. 8, e=0.847 and SW/Sell= 1.489. In Fig. 8, note the uniformity of the intermolecular channel, i.e., the way “bumps pack in hollows”. Since calculations stress non-bonded intermolecular interactions (AEd, Table 3)) it was of interest to compute the non-bonded energy of an atomic probe revolving beyond the generating ellipse. An hydrogen atom was chosen as probe, and Hill-equation parameters were as in MM2 [41]. The chosen trajectory was an ellipse (external ellipse in Figs. 5-7), defined by producing the semiaxes of the generating ellipse by a distance 1.At Z= 0.2 nm, the energy is highly repulsive. At 0.25 nm, it is mildly attractive. In Figs. 5-7, the external ellipse
94
is drawn at 1=0.3 nm, less than the order of intermolecular internuclear distances in the dimers (closest approaches N 0.4 nm). The computed energy on the trajectory is traced as a thick curve [ 421, where the trajectory serves as the coordinate zero. Thus, points within the external ellipse (all points in Figs 57, where 1= 0.3 nm) correspond to attraction. Unlike tracings of the molecular electrostatic potential (MEP), the curve obtained for the non-bonded energy is unstructured: Fig. 5 cannot predict the constellations in Figs. 6 and 7. One cannot claim that, in this context, shape is of minor significance [ 43 1. Rather, unlike in computing MEPs, a probe cannot be representative of an actual environment in a neighboring molecule. On replacing the H-probe by, say, a Cl-atom in a second molecule, two parameters change in the Hill equation (r* and E [8,34] ) and the computed energy is not linear in either. Finally, one is curious to know why electrostatics play such a minor role. To this end we inspected the electrostatic energy of a revolving unit of positive charge. In Fig. 9, this is done for the molecule of Fig. 5, again at a distance of 0.3 nm from the generating ellipse. The charge distribution has been estimated by the modified Del Re method [44]. The ensuing curve does have structure, and minima (arrowheads attracting the probe) are well defined. Maxima (humps repelling the probe) are however weak and shallow. Unlike the nonbonded energy, the electrostatic energy is proportional to the energy of the probe. For an actual atom in a neighboring molecule, charged as it actually is (Qcl,- 0.168; qn, 0.4 - 0.6 electron units), all fine structure would level off. INTRAMOLECULAR
CHANGES
Now to descriptors of the second type. The question to be answered is as follows: when brought to a given environment, does a given molecule contract or expand? Again we consider the van der Waals’ body. Atoms closer to each other than the sum of assigned radii form an overlap lens. Its volume is a measure of their proximity. It is meaningful to distinguish two types of lens: (a) formed by bonded atoms and reflecting bond-lengths; (b ) formed by non-bonded atoms and reflecting bond angles or dihedral angles. Usage of van der Waals’ radii is advantageous, since atoms removed from each other farther than A and D in a bonded sequence A-B-C-D do not overlap. Let On be the sum of bonded overlaps and ONBthe sum of non-bonded overlaps. In analogy with the definition of energy gains, define now “overlap gains”: don=o,-c(i)oO,,~ dONn=ONn-Z(i)O”Nn,i (for bimolecular association, i = 1,2). Positive values express an overall con-
95
TABLE 4 Boeyens’ analysis for compounds of stereochemical interest 0,
W-L
se
0 NB
E
ccl CaH&le, eq ax
39.600 132.499 132.483
9.219 9.142 53.696 53.654
- 640.5 - 639.7 - 1898.8 - 1896.9
Decalin, trans
187.532 187.503
83.691 83.506
- 2615.3 -2612.3
43.208 43.187
10.892 10.895
- 605.6 - 605.0
39.623
CiS
ClCH&H$l,
a iz
Overlaps of van der Waals’ spheres in nm3 x 10e3 (i.e., A 3); energies of Boeyens’ overlaps in kcal mol-‘.
traction of the associated components; negative values express an overall expansion. Atomic overlaps were studied by Boeyens, who also used them to evaluate bond energies [ 45,461. Clearly, the usefulness of the tool goes beyond the few computational examples here at hand. To put numbers into perspective, some results of stereochemical interest are given in Table 4. All refer to computed geometries (alkanes by MM2, chlorides as in ref. 3). The overlaps are by van der Waals radii [lo], and the total energies E by Boeyens radii [ 451. On passing from a less to a more congested stereoisomer, dOs and LION are negative, reflecting expansion. This is because bonds generally stretch and valence angles open up. Actual numbers are small. In the two decalins, for example, d (0s + ONB) is only 0.214 x 10m3nm3 molecule-‘. This is because interatomic distances change only slightly when compared with atomic radii, and the overlap volume depends only weakly on the interatomic distance. In fact, the overlap formula [45] leads to [d(overlap)]/[d(distance)]
=n (; d+r)(+)
where d is the distance, and r is the radius assigned to the overlapping atoms. Now referring to the dimers, in all cases studied (butane, pentane, neopentane, cyclohexane, MeCHClMe, Me3CX, MeCHClCHClMe ) , computed geometries correspond to contraction. It is as if intramolecular squeezing enhances the chance of intermolecular attraction. The numbers, again, are small. For the best au-constellation of (MeCHClCHClMe), (Fig. 6), d(On+ONB)=0.207X10-3nm3molecule-’of the same order as the decalin contraction. The corresponding number for the best gu-constellation, presumably the global minimum (Fig. 7)) is even lower (0.109 x 10S3 nm3 molecule- ’ ) .
dE, also indicates a higher distortion in au (Table 3)) but cannot pinpoint the sense of distortion. CONCLUSION
It has recently been stated that “weak non-covalent bonds involving neutral organic molecules are of vital importance, e.g., in molecular biology, drug design, etc.” [ 471. The statement underlies a curious fashion. Computationalists tend to study molecular association either when the partners are quite complex (e.g., &trypsin with pancreatic trypsin inhibitor [ 21b] ), or when the partners are quite simple (e.g., NH:, with HF, HCI with Nz, NH3 with LiF [48] ). The few examples given above serve to suggest that examination of molecules in between is not unrewarding. In particular, it brings to the fore questions of principle that run the risk of being ignored in dealing with the extremes of the complexity scale. Very small dimers can be treated effectively by quantum chemistry [ 481. Molecular mechanics, though created specifically to meet single molecules, can be applied to more complex dimers. No attempt was made here to optimize the force field, and improvements can be envisaged. As it stands, it did characterize head-to-head orientations of organic halogen compounds, the intramolecular response to intermolecular interactions, and the kinking of carbon-chains in condensed phases. As for descriptors of size and shape, it is true that their content of information becomes useful mainly in the context of complex molecules [491. But it is through application to molecules of lower complexity that they can be tested to advantage [ 41. Here we considered descriptors that reflect the overall change in geometry, eccentricity and surface roughness, on going from a single molecule to an associated pair.
REFERENCES 1 2 3 4 5 6 7 8 9 10 11
A.Y. Meyer, J. Chem. Sot. Perkin Trans., 2 (1986) 1567. A.Y. Meyer, J. Comput. Chem., 9 (1988) 18. A.Y. Meyer, J. Mol. Struct., 94 (1983) 95, and previous papers in the series. A.Y. Meyer, Chem. Sot. Rev., 15 (1986) 449. I. Hargittai and M. Hargittai, in Molecular Structures and Energetics, Vol. 2, VCH, Weinheim, 1987. N. Ramasubbu, R. Parthasarathy and P. Murray-Rust, J. Am. Chem. Sot., 108 (1986) 4308. N.L. Allinger, M.A. Miller, F.A. Van Catledge and J.A. Hirsch, J. Am. Chem. Sot., 89 (1967) 4345. N.L. Allinger,M.T. Tribble, M.A. Mil1erandD.H. Wertz, J.Am. Chem. Soc.,93 (1971) 1637. J. Gajdoli and T. Bleha, Collect. Czech. Chem. Commun., 48 (1983) 71. A. Gavezzotti, J. Am. Chem. Sot., 105 (1983) 5220. J.O. Hirschfelder, J. Chem. Phys., 8 (1940) 431.
97 12 13 14 15
16
17 18 19 20 21
22 23 24 25 26 27
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
J.L. Letelier, Y.-N. Chiu and G. Zon, J. Mol. Struct., 63 (1980) 273. U. Burkert and N.L. Aihnger, J. Comput. Chem., 3 (1982) 40. K.B. Wiberg and R.H. Boyd, J. Am. Chem. Sot., 94 (1972) 8426. Since the terms “conformation” and “configuration” have specific meanings in the context ofsingle molecules, we adopt Prelog’s “con&elk&ion” when referring to the mutual orientation of two molecules. See V. Prelog, J. Chem. Sot., (1950) 420. This exotic energy unit was chosen because enthalpies of vaporization, relevant to the discussion and for comparison, are traditionally tabulated in kcaI mol-‘. Excepting its first occurrence and in tables, the unit is abbreviated as kcai mol-‘. E. Hirota, Y. Endo, S. Saito and J.L. Duncan, J. Mol. Spectrosc., 89 (1981) 285. A. Bondi, J. Phys. Chem., 68 (1964) 441. V. Cody and P. Murray-Rust, J. Mol. Struct., 112 (1984) 189. 0. Bastiansen, K. Kveseth and H. MsllendaI, Topics CUR. Chem., 81 (1979) 99. Complementarity has been studied extensively by G. Naray-Sxab& See, for example, (a) G. Nbray-SxabQ J. Am. Chem. Sot., 106 (1984) 4584; (b) G. Naray-Szab6 and P. Nagy, Enzyme, 36 (1986) 44; (c) P. Nagy, J.G. Angyan and G. Naray-Sxab6, J. Mol. Struct., 149 (1987) 169. M. Santavjr and J. Kypr, J. Mol. Graphics, 2 (1984) 47. A.I. Kitaigorodski, Chem. Sot. Rev., 7 (1978) 133. D. Tabor, Gases, Liquids and Solids, 2nd edn., Cambridge University Press, Cambridge, 1979. (a) p. 153; (b) p. 177. R.C. Weast (editor), Handbook of Chemistry and Physics, 52nd edn., The Chemical Rubber Co., Cleveland, OH, 1971-1972,~. D151. T.L. Hill, Lectures on Matter and Equilibrium, Benjamin, New York, 1966, p. 54. F. Reif, Fundamentals of Statistical and Thermal Physics, McGraw-Hill, New York, 1965, p. 251. I am indebted to Professor S. Baer, of the Hebrew University, for bringing this analysis to my attention. L.R. Pratt, C.S. Hsu and D. Chandler, J. Chem. Phys., 68 (1978) 4202. P.E. Schoen, R.G. Priest, J.P. Sheridan and J.P. Schnur, J. Chem. Phys., 71 (1979) 317. F.H. Westheimer, in MS. Newman (Ed.), Steric Effects in Organic Chemistry, Wiley, New York, 1956, pp. 523-555, in particular p. 550. A.A. Bothner-By and C. Naar-Colin, J. Am. Chem. Sot., 84 (1962) 743. R.J. Abraham and E. Bretschneider, in W.J. Orville-Thomas (Ed.), Internal Rotation in Molecules, Wiley, London, 1974. A.Y. Meyer, J. Comput. Chem., 2 (1981) 384. A.Y. Meyer and F.R.F. Forrest, J. Comput. Chem., 6 (1985) 1. A.Y. Meyer, N.L. Allinger and Y. Yuh, Israel J. Chem., 20 (1980) 57. T. Rydiand and R. Stralevik, J. Mol. Struct., 105 (1983) 157; R. Stelevik, J. Mol. Struct., 109 (1984) 397. J.T. Edward and P.G. Farrell, Can. J. Chem., 53 (1975) 2965. A. Raudino, F. Zuccarello and G. Buemi, J. Chem. Sot., Faraday Trans. 2,77 (1981) 1331. See, for example, W.L. Duren, Calculus and Analytic Geometry. Xerox College publishing, Lexington, MA, 1972, p. 700. J.L. Finney, J. Mol. Biol., 119 (1978) 415. N.L. Allinger, Operating Instructions for MM2 and MMP2 Programs, updated as of August 1985, University of Georgia, 1985. Curve-thickening was achieved by the routines of P. Senn, Comp. Chem., 10 (1986) 219. J.S. Rowlinson, Discuss. Faraday Sot., 49 (1970) 30. A.Y. Meyer, J. Chem. Sot., Perkin Trans., 2 (1982) 1199. J.C.A. Boeyens, J. South Afr. Chem. Inst., 26 (1973) 94. J.C.A. Boeyens, South Afr. J. Chem., 33 (1980) 14, and previous publications in the series.
47 48
49
E. Weber, in the Preface to Top. Curr. Chem., Vol. 140, Springer-Verlag, Berlin, 1987. For very recent examples of ab initio work on dimers, see: G.J.B. Hurst, P.W. Fowler, A.J. Stone and A.D. Buckingham, Int. J. Quantum Chem., 29 (1986) 1223; B. Silvi and M. Allawena, J. Mol. Struct., 135 (1986) 225; Z. Latajka, K. Morokuma, H. Ratajczak and W.J. Orville-Thomas, J. Mol. Struct., 135 (1986) 429. D.E. Walters and A.J. Hopfinger, J. Mol. Struct., 134 (1986) 317.