Faraday Discuss., 1997, 106, 63»77 Computer simulations of organic solids and their liquid-state precursors A. Gavezzotti Dipartimento di Chimica Strutturale e Stereochimica Inorganica, Universitaœ di Milano, via V enezian 21, 20133 Milano, Italy Representative calculations on 2-pyridone in its crystalline, liquid and solvated form, show the capabilities and the shortcomings of computer simulation techniques for the description of organic molecular systems.Both ì static œ and ìdynamicœ computations, meaning without and with the explicit treatment of thermal motion, are employed. Relationships of liquid and solution systems to the structure and thermodynamics of the crystalline solids are investigated as far as possible, with regards to efficient and systematic crystal structure prediction from molecular structure.The thesis to be discussed here is that a reliable theory for the prediction of the solid-state structure of small organic molecules must depend on an accurate thermodynamic and kinetic evaluation of the pathways for evolution from disorder to order, including liquid structuring and nucleation. Dynamic methods, including thermal motion, are indispensable for the accomplishment of this task, which remains today a faraway goal in structural physical chemistry. 1 Introduction When Johann Joachim Winckelmann wrote his ìThoughts on the imitation of the works of the Greeks in painting and sculptureœ in 1775, he inspired the artistic tastes of his times, exalting the uncontaminated beauty of statues and temples as centuries had delivered them to him, with a polished appearance of crystalline whiteness.It was later discovered, that these works of art were originally adorned with appendages of all kinds, in vivid colours, which had faded with time. In a similar fashion, X-ray crystallography has been producing, over the last thirty years, a ìneoclassicœ view of molecular crystals, whereby a static framework structure is seen, as if impenetrable to the disruptive action of thermal motion.An X-ray diÜraction experiment has a timescale of days or, in more recent times, of hours, but always some 1016 times longer than the timescale of molecular motion. Structural chemistry, as it emerged and is still emerging from the use of crystallographic databases, is often oblivious of the massive time averaging that these data incorporate, and relies on sometimes quite evanescent structure»property relationships.Thermal motion is vital in understanding the thermodynamic and kinetic behaviour of chemical systems, but, surprisingly, thermodynamics and kinetics have been put aside more and more by structural chemists, in their dream of reconstructing equilibrium properties and reactivity from idealized atomic positions alone.Nowadays, several analytical techniques are capable of resolution on a molecular scale in the analysis of crystal properties. Important, if humble and sometimes underestimated, experimental information comes from calorimetry : although no detail is available on the nature of the energies and their distribution among relevant thermal modes, the enthalpies of phase change and the speci–c heats of the phases are directly comparable to calculated energies and their derivative.This paper is an essay in the application 6364 Computer simulations of organic solids Table 1 Structural and energetic parameters of 2-pyridone crystals (a) Calculated structures packing space density energy group a/” b/” c/” a/degrees b/degrees c/degrees /g cm~3 /kJ mol~1 P21 5.35 5.77 7.23 » 74 » 1.47 80.3 6.58 3.35 11.44 » 60 » 1.44 79.7 P1- 6.08 3.38 11.50 93 73 77 1.45 80.2 6.16 3.59 11.48 63 90 80 1.42 79.7 6.09 3.66 11.10 105 73 82 1.41 80.0 6.18 3.41 10.78 96 97 77 1.44 80.0 5.98 3.48 11.46 105 106 81 1.43 79.5 P21/c 11.43 3.36 13.18 » 120 » 1.44 79.5 9.94 3.35 13.18 » 95 » 1.44 79.5 11.32 5.88 6.87 » 71 » 1.46 81.2 P212121 13.78 5.37 5.76 » » » 1.48 81.6 18.89 6.63 3.56 » » » 1.42 79.5 Pbca 11.07 5.19 15.73 » » » 1.40 77.6 11.41 11.03 7.26 » » » 1.38 75.6 (b) Experimental structures, P212121 T a/” b/” c/” a/degrees b/degrees c/degrees room 13.645 5.900 5.692 » » » 1.379 77.6 120 K 13.564 5.795 5.604 » » » 1.434 78.2 relaxed 13.82 5.41 5.75 » » » 1.47 81.4 of some computer techniques being developed for the analysis of molecular crystalline properties and, more generally, of the phase behaviour of organic substances.Beginning with crystal packing analysis,1 and an evaluation of packing energies,2 one may proceed to the construction of possible crystal structures starting from molecular structures,3 given suitably calibrated intermolecular potentials4 (Ref. 1»4 provide an ample historical perspective on these methods.) As will be discussed, the signi–cance of computer studies of organic solids is largely enhanced through coupling with simulations of the concomi- Fig. 1 The crystal structure of 2-pyridone7 (space group Oxygen atoms are in black. P212121). Solid lines denote hydrogen bonds (O… … …H distance 1.811 CwO… … …H and O… … …HwN angles 135 ”, and 160°, CwO… … …HwN torsion 109°).A. Gavezzotti 65 tant liquid (either pure or solution) states.Primarily because of this, molecular dynamics is the most appropriate computational tool,5,6 as in principle it should encompass both the thermodynamics and the kinetics of the bulk systems. However, it has its pitfalls : a macroscopic one being that the entropy is not directly accessible.Merits and other shortcomings will also become apparent from our results presented here, no doubt providing a prominent topic for the ensuing Discussion. A suitable candidate molecule to illustrate the capabilities of the above-mentioned techniques is 2-pyridone [2(1H)-pyridinone], as easily parameterizable system, for which a certain amount of structural and thermodynamic information is available.This paper will hopefully enable the reader to evaluate what can nowadays be done, using a computer, to describe and predict the phase behaviour of small organic molecules, and, at the same time, highlight the (many) points that still have to be clari–ed and systematized before these methods can be con–dently incorporated in the physical chemistœs toolbox. 2 2-Pyridone : structural and thermochemical data 2-Pyridone (mp\380 K, bp\553 K) can exist in either the lactam, 1, or lactim, 2, tautomeric forms, the former being found in the solid state7 and favoured by solvation.8 All our calculations and simulations refer to the lactam form, as the activation energy for isomerization is in the 26»52 kJ mol~1 range.8 The heat of combustion has been measured,9 and the ensuing calorimetry showed that the enthalpy diÜerence between the two gas-phase tautomers (2.6 kJ mol~1) is within the experimental uncertainty of the measured sublimation enthalpy (86.6 kJ mol~1).The vibrational frequencies and ideal gas thermodynamic functions of 2-pyridone are also known.10 Scheme 1 The formation of a hydrogen-bonded, cyclic centrosymmetric dimer seems an obvious option for the lactam form, boundary conditions allowing ; dimerization energies have been calculated (by CNDO) as 43.4 kJ mol~1,11 or measured as 35.9 or 65.3 kJ mol~1 in benzene;12 the discrepancy between these last values re—ecting the awkwardness of the measurement. Despite this, the crystal structure of pyridone7 (Table 1 and Fig. 1) adopts the so-called catemer motif, an in–nite chain of hydrogen-bonded molecules in a non-centrosymmetric space group. 3 Static calculations ; packing energies A molecular model for 2-pyridone was taken from the neutron diÜraction study performed at room temperature7 (Cambridge Structural Database ref. code PYRIDO02), renormalizing the NwH distance (experimental 1.03 at 1.00 as required for many ”) ”, subsequent applications : the molecule is —at.The packing energy calculated by standard potential parameters4 is 6% lower than the experimental sublimation enthalpy. The structure was relaxed13 by optimizing cell parameters and molecular orientation under the action of the potentials ; the only signi–cant drift was a 8% shortening in the b cell parameter (Table 1). 4 Static calculations ; crystal structure ìpredictionœ? The PROMET computer package3,14 was used. Brie—y, the procedure consists of assuming a rigid molecular model, building from it a few basic dimeric or oligomeric66 Computer simulations of organic solids aggregates (using inversion centres or screw axes), and then combining these units using the action of other symmetry operators, eventually packing the molecule into the most common space groups for organic compounds: P1-, Pbca.Aggre- P21, P21/c, P212121, gation energies (including hydrogen bonding) are calculated by standard optimized intermolecular potentials4 in the 6-exp form, i.e., without explicit coulombic contributions. The merits and pitfalls of such a parameterization have by now been discussed ad nauseam.1h4 A centrosymmetric dimer was –rst constructed (Fig. 2). Its binding energy, 38 kJ mol~1, falls within the boundaries of experimental values. This dimer was used as a starting building block for P1-crystal structures, it was combined with a screw axis to generate layers and structures, and it was combined with two screw axes to access P21/c space group Pbca.A molecular ribbon along a screw axis was used as the building block for crystal structures, and, after combination with a second screw axis, for P21 P212121 crystal structures. After extensive searches of the crystal energy hypersurface, a few thousand generated crystal structures were grouped and sorted according to space group symmetry and similarity in packing mode and cell dimensions, after Niggli cell reduction. Table 1 collects the most cohesive.The centrosymmetric dimer packs in two diÜerent space groups [Fig. 3(a) and 3(c)] with distortion to a ladder arrangement with diÜerent step heights (from 0 to ca. 2 ”) ; similar arrangements are obtained along a screw axis [Fig. 3(b) and 3(d)]. The second, radically diÜerent packing mode (Fig. 4) stems from the catemer motif, in space groups and the X-ray crystal structure has the packing mode of Fig. 4(b), and that P21 P212121; structure was promptly identi–ed by the search procedure, the eÜectiveness of which is thus again con–rmed. A third packing mode is the benzene-like structure of Fig. 5, in space group Pbca, where the centrosymmetric dimer packs essentially undistorted from the starting structure of Fig. 2. The energy results (Table 1) are, as usual, puzzling. An optimistic viewer might claim that complete success in structure prediction has been achieved, since the X-ray crystal structure is the most cohesive among those identi–ed in the search. However, the packing energies of the generated crystal structures (Fig. 6) form a continuum (as usual15,16) between 70 and 80 kJ mol~1, and polymorphic forms are experimentally known to diÜer in energy by as much as 5»10 kJ mol~1.17 Besides, energy diÜerences, Fig. 2 The planar centrosymmetric dimer generated in the PROMET search (O… … …H distance 1.832 CwO… … …H and O… … …HwN angles 125 and 175°) ”,A. Gavezzotti 67 Fig. 3 Dimer or pseudo-dimer packing motifs in the crystal structures generated during the PROMET search (Table 1).Solid lines : hydrogen bonds, lengths of which are given (below) in parentheses in (a) P1-structure (1.857), dashed line : translation of 3.59 (b) structure ”. ”; P21 (2.044 and 2.067), screw axis running perpendicular to the molecular planes ; dashed line : translation of 3.55 (c) structure (1.873), CwO… … …HwN torsion angle 137°; (d) struc- ”; P21/c P212121 ture (1.873), dashed line : pseudo-hydrogen bond, 2.54 ”.Fig. 4 Catemer motifs for structures generated in the PROMET search (codes as in Fig. 3). (a) P21 structure (1.868) ; (b) structure (1.859), virtually indistinguishable from the X-ray crystal P212121 structure of Fig. 1.68 Computer simulations of organic solids Fig. 5 Pbca crystal structure generated by the PROMET search (Table 1).Centrosymmetric dimers are nearly coplanar with a O… … …H distance of 1.828 ”. even in the sorted structures of Table 1, are in the range of fractions of 1 kJ mol~1, while the energetic resolution of the entire method is far from being that high. A more realistic interpretation is that at least –ve crystal structures in four diÜerent space groups must be considered equally probable on purely energetic grounds.In such a situation, dynamic (thermal and kinetic) factors at the –rst stages of molecular recognition and of Fig. 6 Crystal structures generated in the PROMET search in –ve space groups. For graphic clarity, structures have been considered equal if the energy diÜerence was less than 0.5 kJ mol~1 and the cell volume diÜerence was less than 0.1 irrespective of symmetry.The subpopulation ”3, at high density and low packing energy is composed of (unrealistic) non-hydrogen bonded crystal structures.A. Gavezzotti 69 crystal nucleation must be of paramount importance in making the –nal decision on which structure is most likely eventually to appear, and that is the reason why simulations of liquid precursor phases must be extremely important.For the above reasons, the whole matter is still in epistemological limbo18 (we donœt know what we do know and what we donœt know). Much thought will be needed before crystal structure prediction can be said to be a realistic possibility, or to reconcile the above results with ideas put forward by Desiraju in a recent paper19 on crystal engineering (ì mischief, thou art afootœ). 5 Molecular dynamics simulations The GROMOS package20 has been used to perform energy minimization (EM) and molecular dynamics (MD) calculations on the pyridone dimer in the gas phase and in solution, and on the crystalline and liquid phases. GROMOS non-bonded CCl4 Lennard-Jones functions were used, with and CwH groups in the united atom CCl4 approximation; CxO and NwH fragments were described as one charge group each, with a charge separation of 0.38 and 0.28 e, respectively.Bond distances were constrained to experimental values (SHAKE algorithm). Angle bending force constants were all set equal to 200 kcal mol~1 rad~2, and equilibrium angles were the experimental values. A combination of standard GROMOS proper and improper dihedral potentials was used.No intramolecular non-bonded contributions were included. The intramolecular force –eld is anyway scarcely relevant, and is only a means of preventing distortions from a rigid planar conformation. Periodic boundary conditions were imposed, with constant temperature and pressure restraints (coupling constants of 0.1 and 1 ps, respectively, unless otherwise speci–ed). The cut-oÜ distance in all energy calculations was 10 presumably just enough to ensure proper convergence of electrostatic terms in ”, the charge-group approximation (such a cut-oÜ would be hopelessly too drastic in the point-charge approximation). 5.1 The gas-phase dimer The gas-phase dimer cohesive energy at 1 K (i.e., in the absence of thermal energy) is 38 kJ mol~1, of which 36 kJ mol~1 are from electrostatic and 2 kJ mol~1 from Lennard- Jones contributions. In spite of the widely diÜering force –elds, this result is remarkably similar to that obtained using PROMET calculations.On increasing the temperature of the simulation, the average hydrogen-bond distance increases to 1.96 and 2.04 at 100 ” and 150 K, respectively, and already at 200 K the hydrogen bonds break.Nevertheless, the two molecules do not drift apart, but form a sort of stacked, non H-bonded dimer; in Fig. 7 one sees the O… … …H distances increasing, and the distance between centres-ofmass decreasing towards the minimum separation between stacked planes. This con–guration still has a substantial cohesive energy [(14]8) kJ mol~1, electrostatic and Lennard»Jones, respectively]. These results are only indicative, since equilibrium with the lactim form was not considered. 5.2 The dimer in carbon tetrachloride solution Starting atomic coordinates for 231 solvent and two solute molecules for the simulation of the dimer in were generated by geometrical considerations, and the system was CCl4 subjected to energy minimization –rst.The output of EM was used as the starting point for MD simulations. The coupling constants for the temperature bath were 0.01 ps for the solute and 0.2 ps for the solvent ; in spite of the short coupling constant, the temperature of the solute was always ca. 100 K higher than that of the solvent. At 300 K,70 Computer simulations of organic solids Fig. 7 Molecular dynamics simulation of the gas-phase dimer of 2-pyridone at 200 K: (a) distance between centres-of-mass ; (b), (c) O… … …H distances quick disruption of the dimer is observed, and the two molecules de–nitely drift apart in the solution.To obtain signi–cant information, it was pragmatic to decide to lower the temperature: at 200 K the system stabilizes with an overall density of 1.53 g cm~3 (experimental value for carbon tetrachloride at 298 K of 1.59 g cm~3) and a cohesive energy of the solvent of 19 kJ mol~1 (experimental heat of vaporization of 21 kJ mol~1 at 298 K).The solute»solvent cohesive energy (roughly, the solvation enthalpy) is [57 kJ mol~1; no experimental counterpart is available for this last quantity ; small amides have solvation enthalpies of between [40 and [80 kJ mol~1 in water.21 At 200 K, the dimerization energy is 29 kJ mol~1.The dimerization energy of caprolactam, computed with the OPLS force –eld,22 is 59 kJ mol~1 in chloroform at 300 K; thus, computed values span about the same range as the experimental ones (see Section 2). More relevant to the present purposes is the dynamic behaviour of the dimer in solution. At 200 K, one O… … …H hydrogen bond is always preserved, while the other exhibits ìcatemer jumpsœ (as also observed in the dimer of tetrolic acid23) characterized by a sudden increase of the O… … …H distance to ca. 6 (see Fig. 8), resulting in a destabi- ” lization by ca. 12 kJ mol~1, or about half the total electrostatic dimerization energy. Fig. 8 Molecular dynamics simulation of the 2-pyridone dimer in carbon tetrachloride at 200 K: (a), (b), (c) as in Fig. 7. Small circles denote instantaneous breaking of one O… … …H bond. The dashed peak around 90 ps is a longer catemer jump (see Fig. 9).A. Gavezzotti 71 Over and above the numerous approximations and the theoretical difficulties inherent to the computational methods, this result should be accepted as both reasonable and informative, since these ì—ying catemerœ solution structures (Fig. 9) are likely to be responsible for the seeding of the catemer nuclei that eventually lead to the observed crystal structure of 2-pyridone. 5.3 Molecular dynamics of the pyridone crystal MD simulations for the pyridone crystal were performed using the experimental X-ray crystal structure (PYRIDO02; 128 molecules, 2]4]4 unit cells) as a starting point.A modi–ed24 version of GROMOS, more —exible for crystal calculations, was employed. After preliminary tests, isotropic pressure scaling was settled upon, because a largely negative y pressure component resulted, and anisotropic scaling brought about a disturbingly large contraction of the b cell parameter, similar to what was observed in static minimizations (Section 3).Admittedly, all this casts a shadow on the overall reliability of the potentials, but we justify the choice of isotropic scaling as one of the (many) restrictions and approximations on what, after all, is but a computer simulation. The starting temperature was 200 K, and inputs for other temperatures were then taken from output of a preceding simulation at the nearest temperature.Typical run times were 50»80 ps, and inspection of diagrams for the time evolution of energy and structural parameters ensured that the system had reached a steady state at each temperature. Table 2 collects the main results. Rather than following in detail the evolution of intermolecular interatomic distances, centre-of-mass oscillation was monitored by computing the displacement of each of the Fig. 9 The dimer structure around 90 ps in Fig. 8: one O… … …H distance 2.05 second H-bond ”, broken. The molecular planes are nearly perpendicular to one another. Table 2 Results of molecular dynamics simulations for the 2-pyridone crystala [Etot density [Eel [E(l.j.) *H(subl) T /K /kJ mol~1 a/” b/” c/” /g cm~3 /kJ mol~1 /kJ mol~1 /kJ mol~1 5 72 13.62 5.89 5.68 1.386 20 53 75 120 56 13.66 5.91 5.70 1.374 19 53 75 140 53 13.67 5.91 5.70 1.371 19 52 74 170 49 13.70 5.93 5.71 1.365 18 52 73 200 45 13.71 5.93 5.72 1.359 18 52 73 250 38 13.78 5.95 5.75 1.342 18 51 73 300 31 13.87 6.02 5.78 1.313 17 49 70 350 23 14.00 6.05 5.84 1.278 16 47 67 a Owing to isotropic pressure rescaling the a : b : c ratio is constant.72 Computer simulations of organic solids 128 from its reference position (i.e., that in the starting structure).At 170 K, motion is roughly harmonic and oscillation times range from 0.2 to 0.4 ps, corresponding to frequencies of 160»80 cm~1, very reasonable for intermolecular vibrations. For temperatures up to 300 K, centres-of-mass oscillate about an equilibrium position no more than 0.5 away from the X-ray position, a physiological deviation almost certainly due ” to inaccuracies in the potential functions.More interesting is the fact that at 350 K the motions begin to show a marked anharmonicity (Fig. 10). At still higher temperatures, centre-of-mass displacements show a steady increase as the simulation time proceeds (Fig. 11), yielding a dynamic picture of a disgregation of the crystalline edi–ce that can be assimilated to melting, or even vaporization. The fact that this happens at a temperature close to the real melting temperature must be mere coincidence, since the true timescale of melting is such that it could never be observed in such short simulation times.Fig. 12 shows frames after 50 ps runs, at 350 (with substantial disordering, but a retention of the basic backbone) and at 400 K (with complete loss Fig. 10 Molecular dynamics simulation of 2-pyridone crystal at 350 K; displacements of the centre-of-mass of a representative molecule from the ìequilibriumœ (X-ray) position (a, b, c and d denote the x, y, z components and the modulus of the displacement). A marked anharmonic jump is evident near 0.6 ps simulation time.Fig. 11 Molecular dynamics simulation of 2-pyridone crystal at 400 K: each vertical bar is the range of displacements of the 128 centres-of-mass, the solid line is the average over the range. The moire-like shadow is a graphic artifact.A. Gavezzotti 73 Fig. 12 Final frames of the simulation of the 2-pyridone crystal at (a) 120 K, (b) 350 K, (c) 400 K. The horizontal dimensions are 27.3, 28.0 and 28.5 respectively.”,74 Computer simulations of organic solids Table 3 Results of molecular dynamics simulations for the 2-pyridone liquid [Etot density [Eel [E(l.j.) *H(vap) T /K /kJ mol~1 /g cm~3 /kJ mol~1 /kJ mol~1 /kJ mol~1 320 19 1.145 10 42 56 340 16 1.142 9.6 42 56 370 9.7 1.093 8.6 39 52 400 4.6 1.063 8.0 38 51 of short- and long-distance order).Although a thorough analysis of preferential displacements with respect to the average crystal structure is, at present, too complicated, a detailed study of these motions could no doubt cast some light on the molecular mechanisms of crystal disruption, the timescale problem being circumvented by strong overheating. In any case, Fig. 12(b) should be regarded as the ì realistic œ view of crystal packing, as opposed to the ìneoclassicœ view of Fig. 1. 5.4 Molecular dynamics of the 2-pyridone liquid phase A starting model for liquid-phase simulations was obtained by building a box with 108 molecules with centres-of-mass that lay at face-centred cubic lattice nodes with a spacing Fig. 13 Histograms of O… … …N distances in (a) crystal at 300 K, (b) liquid at 320 K.Circles denote the spread over simulation time, the solid line is the average over the spread.A. Gavezzotti 75 of 8 and whose orientation, described by three rotation angles, was taken at random; ”, no prearranged hydrogen bonding was present. Extensive energy minimization transformed this tentative structure into a substantially cohesive one (57 kJ mol~1).MD was started from this optimized structure, with a few preliminary steps at 100 K to assign velocities followed by warm-up cycles to higher temperatures, as for the crystal. Typical run times were 80»100 ps. Energetic and structural results are collected in Table 3. Fig. 13 shows histograms for N… … …O distances in the liquid and in the solid. The one for the crystal has obvious peaks corresponding to its H-bonding pattern.Not unexpectedly, the liquid curve (even for a formally highly supercooled liquid 60 K below the melting point) shows a large reduction of the population of the 2.8 peak (with a 2 : 1 ” ratio of the O… … …N to the O… … …H peak height) and disappearance of the 3.8 peak, ” corresponding to second-nearest distances in the catemer chain.Overall, the liquid curve is not too diÜerent from the computed radial distribution function for Nmethylacetamide. 21 Along with simple chemical intuition, these results show that the pure liquid, aside from a generic propensity for H-bonding, has no built-in preference for long- or medium-range structuring (of course, any such structuring is strongly disfavoured on entropy grounds). Thus, the liquid is an average bath where precursors to the crystal structure lay buried among diÜerent structures and pure disorder.Things may change around the freezing temperature, but the problems of timescale and of potential function accuracy con–ne the dynamical simulation of structuring and nucleation within a freezing liquid as a task for the future, although signi–cant results have already been obtained on phase transitions and melting in small molecular clusters.25 6 Thermodynamics The MD results in Tables 2 and 3 may be used to calculate thermodynamic parameters of phase transitions of 2-pyridone.In doing so, it should be remembered that the total kinetic energy and the intramolecular potential energy in the simulation are somewhat ill-de–ned, since stretching motions are frozen, and the number of potential-energy terms does not match the number of normal modes; besides, the temperature is evaluated from the kinetic energy in equipartition over all degrees of freedom, an assumption which, if legitimate in the continuum limit inherent to the classical approach of MD, is hardly valid, at least for the bending motions.This difficulty in the de–nition of temperature is quite general, and has to be somehow overcome if MD is to be used routinely for thermodynamic purposes.In deriving *Hs and speci–c heats, the following assumptions are eventually made: (a) E(cohesion)\E(electrostatic)]E(Lennard-Jones)/0.97, where the numerical factor accounts for the cut-oÜ in the summations; (b) E(kin, crystal)\E(kin, liquid)\E(kin, gas) ; this is appropriate for intramolecular contributions, while kinetic-energy diÜerences for external modes should amount to a few RT units, but have not been considered in a –rst approximation; (c) *H(subl)\[E(cohesion, crystal)]RT , *H(vaporization)\ E(cohesion, liquid)]RT , *H(fusion)\E(cohesion, liquid)[E(cohesion, crystal), as per assumption (b) ; (d) PdV being negligible for Cp\dH/dT \dE(total)/dT , condensed phases.The calculated crystal density is only 4% lower than experiment, and the slope of the density curve is well reproduced (Fig. 14). The overall shape of the liquid»solid dilatometric curve is quite satisfactory, and the density diÜerence between crystal and liquid at the melting point is 12%, a typical value for organic substances.The calculated cohesive energy signi–cantly underestimates the enthalpy of sublimation; in this respect, simple packing energy calculations gave a much better result, presumably because of the more speci–c parameterization of the crystal potentials.4 The crystal as computed Cp , from the slope of the E(total)/T curve, is 141 J mol~1 K~1, and is, unfortunately, almost temperature independent.For comparison, can be estimated by adding the ideal gas Cp76 Computer simulations of organic solids Fig. 14 Dilatometric curves for 2-pyridone. Full circles : experimental (X-ray) crystal ; half-–lled circles : calculated, crystal ; open circles : calculated, liquid. The vertical bar denotes the melting temperature. contribution, which averages 78 J mol~1 K~1 between 200 and 300 K,10 and the external contribution calculated from a lattice-dynamical treatment of the crystal,26 or 60 J mol~1 K~1, totalling 138 J mol~1 K~1; the agreement with the molecular dynamics result is purely fortuitous, and must arise from extensive cancellation of opposite errors.No thermodynamic data are available for the pyridone liquid. The calculated vaporization and fusion enthalpies (56 and 14 kJ mol~1, respectively) are of the correct order of magnitude (compare with N-methylisobutyramide and N-methylacetamide, *vapH\ kJ mol~1,27 and 56 kJ mol~1,22 respectively ; N-methylacetamide and acetanilide, 67 and 22 kJ mol~1, respectively28). and density of the liquid at 370 K (184 *fusH\9.7 Cp J mol~1 K~1 and 1.09 g cm~3) compare favourably with experimental data for Nmethylacetamide22 (169 J mol~1 K~1 and 0.894 g cm~3 at the same temperature). 7 Final comment Using only a computer, a completely independent determination of the molar volumes and of phase transition enthalpies for an organic compound, even before synthesizing it, could proceed through the following steps : (1) preparation of a molecular model, using standard geometrical parameters, this is easily accomplished, given the wealth of structural information available nowadays; (2) choice of a suitable force –eld ; (3) generation of crystal structures from the molecular model, and obtention of a satisfactory close packing (even if the predicted crystal structure were wrong, the crystal density and enthalpy of sublimation would be approximately right, see Table 1) ; (4) simulation of the liquid phase, either by MD or Monte Carlo, and calculation of the molar volume and of enthalpies of fusion and vaporization.Step 2 is indeed critical, but progress has been recently made in the development of force –elds. Careful consideration, with avoidance of mindless use of default parameters, can lead to reasonable guesses for most organic molecules.Step 3 suÜers from the drawback that the molecule has to be rigid or semi-rigid ; but an overview of possible polymorphism would be obtained. Steps 3 and 4 are sensitive to parameterization, but alsoA. Gavezzotti 77 experimental uncertainties in thermochemistry are typically a few kJ mol~1. The cost : quality ratio of the results would probably be as good as or better than that for experimental determinations, given the cheapness of computers with respect to chemical synthesis and instrumentation. The quantitative treatment of kinetic energy is still unsatisfactory, so that prediction of speci–c heats is more problematic.Entropy is not accessible, so that a thermodynamic prediction of melting and boiling points is still impossible. Much more problematic is a treatment of the kinetics of phase transition. While a dynamic treatment is here by de–nition indispensable, present theoretical methods are still a long way from being up to the task.Some guidelines have been explored in this paper. Thanks are due to Prof. W. van Gunsteren for permission to use the GROMOS package, and to B. P. van Eijck for useful discussions and advice, in the framework of EC Human Capital and Mobility grant ERB-CHRX-CT94-0469.References 1 A. Gavezzotti and G. Filippini, Acta Crystallogr. Sect. B, 1992, 48, 537. 2 A. Gavezzotti and G. Filippini, in T heoretical Aspects and Computer Modeling of the Molecular Solid State, ed. A. Gavezzotti, Wiley and Sons, Chichester, 1997, ch. 3. 3 A. Gavezzotti, J.Am. Chem. Soc., 1991, 113, 4622. 4 A. Gavezzotti and G. Filippini, J. Phys. Chem., 1994, 98, 4831. 5 W. F. van Gunsteren and H. J. C. Berendsen, Angew. Chem., Int. Ed. Engl., 1990, 29, 992. 6 W. L. Jorgensen, Chemtracts: Org. Chem., 1991, 4, 91. 7 U. Ohms, H. Guth, E. Hellner, H. Dannohl and A. Schweig, Z. Kristallogr., 1984, 169, 185. 8 C. Adamo and F. Lelj, Int. J. Quantum Chem., 1995, 56, 645. 9 S. Suradi, N. El Saiad, G. Pilcher and H. A. Skinner, J. Chem. T hermodyn., 1982, 14, 45. 10 K. C. Medhi, Bull. Chem. Soc. Jpn., 1991, 64, 1944. 11 A. Fujimoto and K. Inuzuka, Bull. Chem. Soc. Jpn., 1990, 63, 2292. 12 S. Mille–ori and A. Mille–ori, Bull. Chem. Soc. Jpn., 1990, 63, 2981. 13 D. E. Williams, PCK83, QCPE Program 548; Quantum Chemistry Program Exchange, Chemistry Department, Indiana University, Bloomington, IN, 1983. 14 PROMET(5), A Program for the Generation of Possible Crystal Structures from the Molecular Structure of Organic Compounds, University of Milano, 1995 (available from the author upon request). 15 A. Gavezzotti, Acta Crystallogr. Sect. B, 1996, 52, 201. 16 B. P. van Eijck, W. T. Mooij and J. Kroon, Acta Crystallogr. Sect. B, 1995, 51, 99. 17 A. Gavezzotti and G. Filippini, J. Am. Chem. Soc., 1995, 117, 12299. 18 A. Gavezzotti, Acc. Chem. Res., 1994, 27, 309. 19 G. R. Desiraju, Angew. Chem., Int. Ed. Engl., 1995, 34, 2311. 20 W. F. van Gunsteren, GROMOS, Groningen Molecular Simulation Program Package, University of Groningen, 1987. 21 W. L. Jorgensen and C. J. Swenson, J. Am. Chem. Soc., 1985, 107, 1489. 22 W. L. Jorgensen and C. J. Swenson, J. Am. Chem. Soc., 1985, 107, 569. 23 A. Gavezzotti, G. Filippini, J. Kroon, B. P. van Eijck and P. Klewinghaus, Eur. J. Chem., 1997, 3, 893. 24 B. P. van Eijck, personal communication. 25 L. S. Bartell, J. Phys. Chem., 1995, 99, 1080. 26 G. Filippini, personal communication. 27 J. S. Chickos, D. G. Hesse and J. F. Liebman, J. Org. Chem., 1989, 54, 5250. 28 J. S. Chickos, C. M. Braton, D. G. Hesse and J. F. Liebman, J. Org. Chem., 1991, 56, 927. Paper 7/01436H; Received 28th February, 1997