J. Phys. Chem. B 2001, 105, 10573-10578
10573
Crystal Structure Prediction for Six Monosaccharides Revisited Bouke P. van Eijck,*,† Wijnand T. M. Mooij,†,‡ and Jan Kroon†,§ Department of Crystal and Structural Chemistry, BijVoet Center for Biomolecular Research, Utrecht UniVersity, Padualaan 8, 3584 CH Utrecht, The Netherlands, and Accelrys Ltd., 230/250 The Quorum, Barnwell Road, Cambridge CB5 8RE, England ReceiVed: June 20, 2001; In Final Form: August 20, 2001
Six monosaccharides were studied as prototypes for flexible molecules, forming crystal structures with complex hydrogen bond patterns. In earlier work the prediction of these crystal structures had been attempted using empirical force fields. Now ab initio energies have been calculated for 20 hypothetical structures of each substance, along with corrections to the free energy at room temperature. Five experimental structures corresponded to the global free energy minimum, the sixth ranking second at a relative free energy of only 1.3 kJ/mol. This strongly suggests that a thermodynamic approach to crystal structure prediction will often be adequate, kinetic effects being dominant only in a minority of cases.
Introduction Most crystal structure predictions start with generating a list of hypothetical structures with low energy.1 Then a major problem becomes manifest: many structures are encountered within a limited energy range. Usually the observed polymorphs are present in the list, but not with the lowest energy. This may be due to (1) inaccuracies in the energy calculation, (2) thermal effects, (3) kinetic effects. In recent work we have studied the first two effects for the hydrogen-bonded structures of the two flexible molecules glycol and glycerol. After a preliminary study,2 the energy calculation was carried out as accurately as possible, using quantummechanical methods. The energies of the experimental structures were within 1.5 kJ/mol of the global minimum.3 Next, the effects of harmonic lattice vibrations were found by standard latticedynamical calculations, which allowed a ranking of the structures according to their free energies. In this list glycol ranked first and glycerol second, only 1.1 kJ/mol above the global free energy minimum.4 This is well within the estimated accuracy of the calculations. Thus it seems that kinetics is not a primary factor for the prediction of these structures, crystallized at low temperatures from the liquid phase. This cannot be true in general, since it is known that the choice of crystallization conditions may determine which polymorph crystallizes preferably. So, in fact, observed structures may well be metastable. Yet we thought it worthwhile to extend our studies to larger molecules, crystallized from solutions at room temperature, to see how far the purely thermodynamic approach can be pushed. The substances chosen were the six hexapyranoses for which extensive lists of hypothetical structures were already available from our earlier work, using empirical force fields.5,6 These molecules contain a considerable amount of conformational flexibility. Moreover, the orientation of the hydroxyl groups is mainly determined by the crystal packing, so the conformational * Corresponding author. E-mail:
[email protected]. † Utrecht University. ‡ Accelrys Ltd. E-mail:
[email protected]. § Deceased.
problem cannot be solved independently from the packing problem. In the present work we first determined the minima in the total crystal energy that correspond to the experimental structures and established the correspondence between observed and calculated geometrical features. Second, we compared the crystal energy of these structures to that of other, hypothetical, crystal packings. Ideally, the energy minima of the experimental structures should correspond to the global (free) energy minima. Although an exhaustive study of possible minima was beyond our computational resources, it was hoped that a carefully selected set of 20 structures for each compound would form a representative sample of the low-energy region of the crystal energy hypersurface. Energy Calculation The intermolecular energies were calculated using a potential that was developed earlier from high-level ab initio calculations on methanol dimers and trimers.7 It involves atomic multipole moments, atomic dipole polarizabilities, a damped atom-atom r-6 dispersion contribution, and an exponential repulsion term that is anisotropic for oxygen. The atomic multipoles are derived from an ab initio wave function at the SCF/DZ(2dO) level. The incorporation of polarizability is essential to allow the transferability of the potential from the gas phase to the crystalline state. This intermolecular potential was shown to give satisfactory crystal structure predictions for methanol, ethanol, dioxane, and propane.8 To obtain the total crystal energy, the intramolecular energy should be added. For flexible molecules this is not a trivial problem. In previous work we applied the empirical force field MM3 to find the intramolecular contribution, which gave reasonable total energies for the crystal structures of glycol and glycerol.2 However, in preliminary work on hexapyranoses9 this approach appeared to fail. Both MM3 and OPLS were tried, and it was found in both cases that at least one experimental crystal structure had an unbelievably large relative total energy of over 20 kJ/mol. Moreover, the two sets of energies were quite inconsistent. Therefore, we resorted to a crystal optimization procedure that combines an intermolecular potential with explicit ab initio
10.1021/jp012366j CCC: $20.00 © 2001 American Chemical Society Published on Web 09/29/2001
10574 J. Phys. Chem. B, Vol. 105, No. 43, 2001 calculations for the intramolecular energy.3 In this work our ab initio intermolecular potential7 was used together with an SCF/ 6-31G* treatment of the intramolecular part. As no empirical energy calculations are involved, the thus calculated crystal energies will be referred to as ab initio energies. In short, in the crystal optimization procedure an ab initio energy and first and second derivatives are calculated for the free molecule in the geometry imposed by the crystal surroundings. This yields, within the harmonic approximation, a functional form for the ab initio intramolecular energy and forces around the starting geometry. The total crystal energy and forces are obtained by adding the intermolecular energy and forces as calculated from the intermolecular potential, and this crystal energy is then minimized. In this work minimization was deemed to be complete when the residual forces were smaller than 0.0004 kJ mol-1 Å-1. In the energy minimizations the ab initio calculations have to be redone when intramolecular structural changes become too large, because the description of the intramolecular energy is only valid around the starting geometry in which the ab initio energy and first and second derivatives are calculated. Criteria of 0.1 Å for bond lengths and 5° for bond angles and dihedral angles were chosen as the limits for the applicability of the harmonic approximation. A problem with this iterative method is that sometimes unexpectedly large oscillational behavior occurs and exact convergence may be hard to reach. The noise level was in the order of 1 kJ/mol; for the most promising structures a second cycle of calculation was done. The procedure, which led to successful crystal structure predictions for glycol and glycerol, has recently been described in detail.3 It has been implemented in a local version of the TINKER package,10 which calls GAMESS-UK11 for the ab initio calculations and MOLDEN12 to fit the atomic multipole moments to the electrostatic potential around the molecule. A molecular cutoff radius of 20 Å, based on the centers of geometry, was used throughout all calculations. In the present work we applied this computational strategy to a limited number of hypothetical structures for the six hexapyranoses. Here a new problem may arise: some structures contain intramolecular and intermolecular hydrogen bonds, and it is not obvious whether the cooperativity is adequately modeled in such cases. In the Appendix we estimate that the errors introduced here are limited to 1 kJ/mol for a system of two hydrogen bonds. Selection of Structures Our first study5 involved the standard GROMOS force field13 and Kouwijzer’s all-atom extension for carbohydrates.14 The latter model produced apparently good results for four of the six hexapyranoses. However, in both sets of calculations the experimental structures were found to have an unreasonably high relative energy for β-D-galactose and β-D-glucose. In a more recent study,6 slight modifications to GROMOS resulted in the united-atom force field UNITAT.6 Its performance was compared to the OPLS all-atom carbohydrate force field.15 Table 1 summarizes all previous results; ∆E represents the energy of the experimental structure with respect to the global minimum and R its ranking. It appears that these values are rather sensitive to the force field used. The root-mean-square difference between the energies in UNITAT and OPLS is more than 10 kJ/mol. Thus the rankings in Table 1, which are based on energy differences of that order of magnitude, can hardly be considered as significant. For the present work a list of hypothetical OPLS structures was taken as the basis. In the previous work, these structures
van Eijck et al. TABLE 1: Results of Previous Crystal Structure Predictiona GROMOS5 Kouwijzer5 UNITAT6 OPLS6 substance R-D-galactose R-D-glucose R-D-talose β-D-allose β-D-galactose β-D-glucose
refcode
R
ADGALA 8 GLUCSA 8 ADTALO 20 COKBIN 6 BDGLOS 105 GLUCSA 213
∆E
R
10.5 6.3 10.0 2.9 20.5 23.8
1 3 1 2 88 384
∆E 4.6 3.8 20.5 21.3
R
∆E
R
∆E
7 13 15 28 67 48
5.6 5.5 1.7 9.2 13.0 10.2
9 45 68 65 66 5
8.5 14.1 16.5 12.7 14.7 4.9
a ∆E is the energy (kJ/mol) of the experimental structure with respect to the structure of lowest energy, R is the corresponding ranking.
had been generated with united-atom force fields and subsequently studied in OPLS. Occasionally this approach will miss structures that correspond to high-energy minima in the preliminary force fields, or perhaps not to minima at all. To investigate the importance of this effect, the structure generation was now repeated directly in the OPLS force field, using starting structures where all crystallographic parameters took random values.16 Moreover, in each starting structure random values were also assigned to the O-C-C-O dihedral angle of the exocyclic hydroxymethyl group and to the five hydroxyl C-CO-H dihedral angles. No lower global minima were found, but new structures with intermediate energy did turn up. Therefore, the values of ∆E in Table 1 remained unchanged, but the rankings R deteriorated to values between 6 and 131. Such high rankings indicate that, when these empirical force fields are used for the screening of feasible crystal packings, one should study for each compound at least the 200 most favorable crystal structures to be on the safe side for a fully reliable structure prediction. Unfortunately, our method of ab initio calculations is quite time-consuming: on the average, a full energy minimization for one structure took between one and 2 days on a personal computer equipped with a 550 MHz Pentium III processor. For the present study we chose the most favorable structures in the OPLS force field, together with structures derived from the six most favorable ones in UNITAT. The experimental structures were added, augmented with the structures from the list that came closest to them. Furthermore, care was taken to include structures with each of the three possible conformations for the exocyclic CH2OH group and to include at least one structure with a hydrogen bond that was exclusively intramolecular. All in all, twenty structures, in space groups P21 or P212121, were investigated for each of the six compounds. Obviously, this selection is too limited to ensure the presence of the global minimum or to establish a reliable energy ranking, but it should be possible to obtain an impression of the ∆E values that can be expected from the ab initio force field. Results An important requirement for a force field is that it should reproduce the observed crystal geometries without excessive changes in geometry. Energy minimizations were started from the experimental structures as well as from the corresponding OPLS-optimized structures produced by the crystal structure prediction. For five substances the results (denoted by E and P, respectively) were essentially the same, with almost identical geometries and energy differences within 0.7 kJ/mol. This is just the noise level inherent to the energy minimization algorithm. The problematic case was β-D-allose, where the P-structure was 4.5 kJ/mol lower in energy than the E-structure. The difference is in the intermolecular hydrogen bonding of the hydroxyl hydrogen atom H6. The hydrogen bonding scheme
Crystal Structure Prediction
J. Phys. Chem. B, Vol. 105, No. 43, 2001 10575
TABLE 2: Optimized Crystal Structures for Six Hexapyranosesa ref ADGALA0314 GLUCSA1017 ADTALO1024 COKBIN25 BDGLOS0126 GLUCSA0414
obs calc obs calc obs calc obs calc(E) calc(P) obs calc obs calc
a
b
c
5.90 5.85 4.97 4.85 7.63 7.51 4.92 4.85 4.73 7.70 7.58 6.60 6.42
7.84 7.65 10.37 10.25 8.08 8.16 11.92 11.64 11.69 7.77 7.92 9.01 8.81
15.68 15.66 14.85 14.79 12.10 11.73 12.80 12.88 13.04 12.66 11.88 12.72 12.79
∆X
F
∆θ
0.067
0.8
0.055
1.3
0.051
2.7
0.13 0.24
1.6 3.6
0.019
3.6
0.067
2.7
1.65 1.71 1.56 1.63 1.60 1.66 1.59 1.65 1.65 1.58 1.68 1.58 1.65
∆τCO
∆τOH
5.8
12.0
2.3
6.5
4.8
11.4
3.7 3.9
27.7 13.9
3.3
11.9
5.7
16.2
a
All structures crystallize in space group P212121. Units are Å and degrees. ∆X is the net translation of the center of mass, calculated via fractional coordinates. ∆θ is the net rotation of the molecule and F is the density (g cm-3). ∆τCO is the deviation in the exocyclic O-C-C-O dihedral angle, ∆τOH is the largest deviation in the five C-C-O-H dihedral angles. For β-D-allose, E and P denote starting geometries from the experimental geometry and the OPLS-predicted structure, respectively.
is O6-H6 ‚‚‚O1 in the P-structure and O6-H6 ‚‚‚O5 in the E-structure, where the interaction with O1 is just the minor component of a bifurcated hydrogen bond. In view of the experimental uncertainty in the hydrogen atom position, it was difficult to decide which structure is the correct one, and we opted tentatively for the energetically more favorable P-structure. The geometry shifts that occurred upon the energy minimizations are given in Table 2. The root-mean-square deviation in the cell axes was 2.4%; the calculated densities were between 3.6% and 6.2% larger than the observed ones. Qualitatively this is as expected: the calculated densities refer to static structures, whereas thermal expansion is present in experimental structures even at very low temperatures. The exocyclic O-C-C-O dihedral angles were reproduced quite well; the correspondence for the C-C-O-H hydroxyl dihedral angles was less satisfactory. However, the uncertainties in the experimental hydrogen coordinates are certainly not negligible14 except in the case of R-D-glucose which was studied by neutron diffraction.17 It is encouraging to see that the deviation in the hydroxyl torsional angles (∆τOH in Table 2) is lowest for this compound. The problems with the reproduction of the structure of β-D-allose are reflected in exceptionally large deviations in the center of mass. A very discriminating test of the quality of a force field is that it should be able to predict the observable polymorph(s) of a given substance. The usual criterion for this purpose is simply the energy (E), although the free energy (A) would theoretically be preferable. We have recently developed a program that performs standard lattice-dynamical calculations to obtain the difference between these two quantities, assuming harmonic lattice vibrations.4 These corrections were calculated in the OPLS force field at 300 K. Although R-D-galactose and β-Dglucose were studied at 95 K, the crystals had been obtained at room temperature. Even if another modification would be thermodynamically more stable at liquid nitrogen temperature, it is doubtful whether the necessary phase transformation could actually take place. Thus it seems more appropriate to calculate the free energy at the temperature of crystal growth. The rankings and relative (free) energies of the energyminimized P-structures among the 20 chosen hypothetical structures are given in Table 3. This table also includes a comparison of the ab initio energies with those calculated by the two empirical force fields. It should be emphasized again that for these flexible molecules only the total energy is a useful criterion: structures that are close together in total energy may differ by as much as 40 kJ/mol in their intramolecular energies.
TABLE 3: Ab Initio SCF Relative Energiesa substance
∆E
RE
∆A
RA
D(UNITAT)
D (OPLS)
R-D-galactose R-D-glucose R-D-talose β-D-allose β-D-galactose β-D-glucose
(16.8) (4.8) 3.2 4.5 6.4 1.9
1 1 2 3 3 3
(11.7) (7.3) 2.4 0.5 1.0 (3.6)
1 1 2 2 2 1
9.9 9.3 8.3 6.4 8.0 9.0
12.9 13.7 11.2 10.3 11.0 8.4
a (Free) energy differences ∆E, ∆A (kJ/mol) and rankings R , R E A refer to the experimental structure with respect to the global (free) energy minimum at 300 K, as found from 20 selected structures for each substance. In case of a ranking R ) 1 the entry in parentheses refers to the (free) energy difference with the second best structure. D is the root-mean-square difference between the 20 ab initio energies and the corresponding empirical values.
Obviously, it is not viable to simply limit crystal structure generations to conformations with low gas-phase energies, as intermolecular hydrogen bonds are essential for the understanding of carbohydrate crystal structures. Additional Calculations All experimental structures were found within 3 kJ/mol of the global free energy minima. It is somewhat surprising that such a good result can be obtained at the SCF level, where dispersion is lacking and significant basis set superposition errors are present. Barrows et al.18 have observed that the SCF/6-31G* level produced relative energies for 11 conformers of glucose that were in surprisingly good agreement with the best composite level (mean absolute error 0.8 kJ/mol). However, such studies involve gas-phase minima with the maximum number of intramolecular hydrogen bonds, and the conclusions do not necessarily apply to crystal structures where the variety of hydrogen bonds is larger. Hopefully, here too the errors will cancel to a large extent when the energies of various structures are compared. The largest effect is expected when different conformations of the hydroxymethyl group are involved. Indeed, in the energetically most favorable structures for R-D-talose and β-D-allose this group has a conformation that differs from the experimental one. To investigate this, calculations at higher levels of theory have to be carried out. For glycol and glycerol MP2/6-311+G(2d,2p) calculations at fixed geometry could be afforded.3 Hexapyranoses are, however, too large for such calculations. Therefore, we resorted to DFT calculations, which provide a computationally tractable alternative. Although it is questionable whether dispersion interactions are properly described in density func-
10576 J. Phys. Chem. B, Vol. 105, No. 43, 2001
van Eijck et al.
TABLE 4: Ab Initio DFT Relative Energiesa substance
∆E
RE
∆A
RA
D
N5
R-D-galactose R-D-glucose R-D-talose β-D-allose β-D-galactose β-D-glucose
(14.5) (4.7) (0.7) 2.9 4.7 6.2
1 1 1 5 3 5
(10.9) (7.4) (1.5) (0.4) (0.7) 1.3
1 1 1 1 1 2
4.0 4.0 2.4 3.1 3.5 3.1
1 1 2 5 3 2
a See Table 3 for explanation. D is now the root-mean-square difference between the intramolecular energies from SCF/6-31G* versus PW91/EXT calculations. N5 is the number of structures within a free energy window of 5 kJ/mol.
tional theory, it has been shown that the exchange and correlation functional PW91 of Perdew and Wang19 at least partially recovers the attraction in van der Waals bonded systems,20 in contrast to the BLYP and B3LYP functionals. Calculations of the intramolecular energy were performed without further geometry optimization, using DMol3 (version 4.2) within the Cerius2 suite of programs.21 DMol uses numerical basis functions that are obtained from solution of the DFT equations for the individual atoms.22 Such basis sets are very accurate and minimize superposition effects; in this work the so-called extended basis set (EXT) was used, a doubled numerical basis set with one set of p- and two sets of d-polarization functions. The results (Table 4) are quite promising. Except for β-D-glucose, the relative energies improved to the extent that the observed structures corresponded to the global free energy minima. Taking the PW91/EXT level as reference, Table 4 shows that the SCF/6-31G* energies exhibit a rootmean-square deviation of a few kJ/mol. Discussion It is seen in Table 2 that the ab initio force field reproduces most features of the experimental geometries satisfactorily. Considering that the intermolecular potential was derived from free dimers and trimers of methanol, this is a nontrivial accomplishment. Empirical force fields do not perform significantly better. Recently Williams23 concluded that no static model can be expected to reproduce the observed cell edges to greater accuracy than about 3%. Moreover, in his work each intramolecular geometry was kept fixed. For carbohydrates, the adoption of the observed hydroxyl group orientations will virtually fix the crystal structure. As regards structure prediction, Tables 3 and 4 suggest that the free energy is a better criterion for crystal structure prediction than is the static energy. Of course, this is only as it should be, but the calculation of the difference A - E may not be very accurate: it is done using an empirical force field (OPLS) and neglecting all effects of anharmonicity. Nevertheless, the limited amount of experience obtained so far indicates that it is worthwhile to apply the correction. The rankings presented in Tables 3 and 4 have to be viewed with care, since the number of structures was limited to 20 for each compound. With that proviso, they appear to be entirely satisfactory. If we estimate the accuracy of the free energy differences to be about 5 kJ/mol, a small number of structures could correspond to the global minimum. This number is given as N5 in Table 4. For R-D-galactose and R-D-glucose, the second best structure is so high in free energy that the difference can be considered as significant. In the other four cases the experimental structure ranks either first or second, with a small value for ∆A: under 2.4 kJ/mol in the SCF approach (where the structures were consistently optimized) and 1.3 kJ/mol in the subsequent DFT calculations (where the minimized geom-
etries were taken over). All unreasonably high relative energies obtained previously (∆E in Table 1) have disappeared. If kinetic effects would be prominent in the crystallization of one of these six carbohydrates, we might have encountered an unacceptably large value of ∆A. This was not the case, which demonstrates that calculation of the thermodynamically most stable polymorph is a useful approach to crystal structure prediction. A comparison of ab initio energies with the values from OPLS and from UNITAT (Table 3) shows large discrepancies. The root-mean-square differences for intra- and intermolecular energies separately are also in the order of 10 kJ/mol. Surprisingly, the simple united-atom force field UNITAT reproduces the ab initio energies somewhat better, with β-D-glucose as the exception. A considerably better agreement between OPLS and 6-31G* (rms difference 3 kJ/mol) is found from the energies reported by Barrows et. al18 for gas-phase conformations. In crystal structures, the balance between inter- and intramolecular forces is delicate and the extent to which energetically unfavorable conformational changes are compensated by a better packing energy is dependent on the force field. This means, unfortunately, that the preliminary ranking obtained with empirical force fields is not at all to be trusted. It is even quite conceivable that for one or more molecules the global minimum has been missed in this work. If the present exercise had been carried out without knowledge of experimental information, four out of the six observed structures would not have been included in a test set of 20 structures. Thus no genuine structure predictions could be made, since the number of structures selected for detailed study had to be limited. As the availability of cheap computer power still appears to increase, this may not be a permanent problem. At least as important will be the development of better simple force fields for the generation of possible structures: one needs a small starting set, with reasonable confidence that it contains the experimentally observable polymorph(s), before one can employ computationally expensive methods to evaluate accurate free energy differences. Acknowledgment. We are grateful to Prof. F. B. van Duijneveldt for helpful discussions. This work was supported by the Council for Chemical Sciences of The Netherlands Organization for Scientific Research (CW-NWO), in the framework of the PPM/CMS crystallization project. Appendix: Nonadditivity in Intramolecular Hydrogen Bonds In our modeling approach we use completely independent models for intramolecular and intermolecular interactions. The latter part contains contributions from the direct electrostatic interaction and from first-order intermolecular polarization. For both purposes atomic multipoles are needed, which are derived from the wave function calculated for the isolated molecule. This simple model was seen to provide an excellent description of the nonadditivity in methanol trimers.7 In crystal structures of carbohydrates, chains of hydrogen bonds occur that can be either intra- or intermolecular. Such a situation can be modeled by a chain of methanol monomers. A case of exclusively intermolecular hydrogen bonds is schematically depicted in Figure 1a; the polarizing field at B is the sum of the electrostatic fields of A and C. However, a hydrogenbonded chain that contains an intramolecular hydrogen bond is treated differently, as depicted in Figure 1b where the A‚‚‚B hydrogen bond is intramolecular. Now the polarizing field at B does not contain the field due to A, as the polarization model
Crystal Structure Prediction
J. Phys. Chem. B, Vol. 105, No. 43, 2001 10577
Figure 1. Chain of hydrogen bonds considered as two purely intermolecular ones (a) or as a combination of an intramolecular and an intermolecular hydrogen bond (b).
TABLE A1: Electrostatic and Polarization Interaction Energies (kJ/mol) in Methanol Trimers Calculated as Trimer ABC and as Dimer of AB with Ca trimer
∆Eelec AB -30.97 -18.64 -30.97 -12.98
G A E D
elec pol elec ∆Epol AB ∆E(AB)C ∆E(AB)C ∆EABC
-6.71 -36.65 -2.44 -38.68 -6.71 -23.54 -2.86 -7.97
-8.75 -6.55 -5.44 -3.07
pol ∆EABC
∆∆E
-64.57 -17.28 1.23 -55.34 -11.08 -0.11 -55.85 -11.42 -0.61 -20.38 -5.65 0.85
a
G, A, E and D denote trimer geometries as given in Figure 2. The C monomer is the left molecule in G and E, the middle one in A, and the right/top one in D.
does not contain intramolecular contributions. The influence of the intramolecular hydrogen bond is more indirect: it is present in the calculation of the charge distribution of the AB molecule and is therefore reflected in the atomic multipole moments. In this example, the dipole moments of the A and B hydroxyl groups are enhanced compared to an isolated hydroxyl group, which obviously influences the calculated interaction with molecule C. To assess the errors associated with this approach, we performed some calculations on methanol trimers, in which we considered one of the hydrogen bonds as “intramolecular”. The interaction energy for a trimer is defined as
∆EABC ) EABC - EA - EB - EC We now take molecules A and B as forming a single molecule. We calculate the charge distribution for the AB supermolecule and let AB then interact with C, resulting in ∆E(AB)C. The trimer energy ∆EABC is then calculated as the sum of ∆E(AB)C and the “intramolecular” energy of the AB supermolecule, which is just the interaction energy within the AB dimer (∆EAB). This latter energy should be calculated with the normal, monomer charge distributions for A and B. In our model we assume that both approaches yield the same interaction energies. So, the difference ∆∆E, defined as
∆∆E ) ∆EABC - (∆EAB + ∆E(AB)C) should be zero. Note, however, that the partitioning of the electrostatic and the polarization energy is different (∆Eelec ABC * elec pol pol pol + ∆E , and ∆E * ∆E + ∆E ). Geometries ∆Eelec AB (AB)C ABC AB (AB)C representing four types of methanol trimers were taken from the ab initio work involved in the parameterization of our intermolecular potential.7 We calculated atomic multipole moments for the methanol monomer and the AB supermolecule at the ab initio level employed for trimers in that work (SCF using the IOM basis without f and g functions).
Figure 2. Overview of the methanol trimer geometries that were used for testing the nonadditivity in intramolecular hydrogen bonds. The naming is the same as in ref 7, where all structural details have been deposited. Used here are G: ROO ) 2.85 Å; A: RHO ) 2.20 Å; D: RHO ) 2.20 Å.
Results for the trimer calculations are given in Table A1. It appears that the electrostatic + polarization energy in the (AB)C “dimer” approach indeed gives a reasonable description of the ABC trimer interaction energy. The difference in intermolecular polarization energy is partly compensated for by the difference in electrostatic energy. The remaining values for ∆∆E are in the order of 1 kJ/mol. So, this is the kind of error that can be expected in comparing intramolecular and intermolecular hydrogen-bonded chains. References and Notes (1) Verwer, P.; Leusen, F. J. J. in ReViews in Computational Chemistry; Lipkowitz, K. B.; Boyd, D. B., Eds., volume 12; Wiley-VCH: New York, 1998; Chapter 7. (2) Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Am. Chem. Soc. 2000, 122, 3500. (3) van Eijck, B. P.; Mooij, W. T. M.; Kroon, J. J. Comput. Chem. 2001, 22, 805. (4) van Eijck, B. P. J. Comput. Chem. 2001, 22, 816. (5) van Eijck, B. P.; Mooij, W. T. M.; Kroon, J. Acta Crystallogr. 1995, B51, 99. (6) van Eijck, B. P.; Kroon, J. J. Comput. Chem. 1999, 20, 799. (7) Mooij, W. T. M.; van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G. C. M.; van Eijck, B. P. J. Phys. Chem. A 1999, 103, 9872. (8) Mooij, W. T. M.; van Eijck, B. P.; Kroon, J. J. Phys. Chem. A 1999, 103, 9883. (9) Mooij, W. T. M. Ab Initio Prediction of Crystal Structures. Ph.D. Thesis, Utrecht University, The Netherlands, 2000. (10) Ponder, J. W. TINKER: Software Tools for Molecular Design, version 3.6.; Washington University of Medicine, 1998. TINKER is available via the Internet at http://dasher.wustl.edu/tinker/. (11) GAMESS-UK is a package of ab initio programs written by M. F. Guest, J. H. van Lenthe, J. Kendrick, K. Schoffel, P. Sherwood, and R. J. Harrison, with contributions from R. D. Amos, R. J. Buenker, M. Dupuis, N. C. Handy, I. H. Hillier, P. J. Knowles, V. Bonacic-Koutecky, W. von Niessen, V. R. Saunders and A. J. Stone. The package is derived from the original GAMESS code due to M. Dupuis, D. Spangler, and J. Wendoloski, NRCC Software Catalog, Vol. 1, Program No. QG01 (GAMESS), 1980. (12) Schaftenaar, G.; Noordik, J. H. J. Comput.-Aided Mol. Design, 2000, 14, 123. MOLDEN is available via the Internet at http://www.cmbi.kun.nl/ ∼schaft/molden/molden.html. (13) van Gunsteren, W. F.; Berendsen, H. J. C. GROMOS, Groningen molecular simulation package; University of Groningen, 1987. (14) Kouwijzer, M. L. C. E.; van Eijck, B. P.; Kooijman, H.; Kroon, J. Acta Crystallogr. 1995, B51, 209. (15) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. J. Comput. Chem. 1997, 18, 1955. (16) van Eijck, B. P.; Kroon, J. Acta Crystallogr. 2000, B56, 535. (17) Brown, G. M.; Levy, H. A. Acta Crystallogr. 1979, B35, 656.
10578 J. Phys. Chem. B, Vol. 105, No. 43, 2001 (18) Barrows, S. E.; Storer, J. W.; Cremer, C. J.; French, A. D.; Truhlar, D. G. J. Comput. Chem. 1998, 19, 1111. (19) Perdew, J. P.; Wang, Y. Phys. ReV. 1992, B45, 13244. (20) Tsuzuki, S.; Lu¨thi, H. P. J. Chem. Phys. 2001, 114, 3949. (21) Cerius2 and DMol3 are products of Accelrys Incorporated, 9685 Scranton Road, San Diego, CA 92121-3752. (22) Delley, B. J. Chem. Phys. 1990, 92, 508.
van Eijck et al. (23) Williams, D. E. J. Comput. Chem. 2001, 22, 1. (24) Ohanessian, J.; Avenel, D.; Kanters, J. A.; Smits, D. Acta Crystallogr. 1977, B33, 1063. (25) Kroon-Batenburg, L. M. J.; van der Sluis, P.; Kanters, J. A. Acta Crystallogr. 1984, C40, 1863. (26) Longchambon, F.; Ohannessian, J.; Avenel, D.; Neuman, A. Acta Crystallogr. 1975, B31, 2623.
