The Thermal Stability of Lattice-Energy Minima of 5-Fluorouracil

4298

J. Phys. Chem. B 2008, 112, 4298-4308

The Thermal Stability of Lattice-Energy Minima of 5-Fluorouracil: Metadynamics as an Aid to Polymorph Prediction Panagiotis G. Karamertzanis,† Paolo Raiteri,‡ Michele Parrinello,‡ Maurice Leslie,£ and Sarah L. Price*,† Department of Chemistry, UniVersity College London, 20 Gordon Street, London WC1H 0AJ, U.K., Computational Science, Department of Chemistry and Applied Biosciences, ETH, Zurich, USI Campus, Via Buffi 13, CH-6900 Lugano, Switzerland, and STFC Daresbury Laboratory, Daresbury, Warrington WA4 4AD, U.K. ReceiVed: October 5, 2007; In Final Form: January 7, 2008

This paper reports a novel methodology for the free-energy minimization of crystal structures exhibiting strong, anisotropic interactions due to hydrogen bonding. The geometry of the thermally expanded cell was calculated by exploiting the dependence of the free-energy derivatives with respect to cell lengths and angles on the average pressure tensor computed in short molecular dynamics simulations. All dynamic simulations were performed with an elaborate anisotropic potential based on a distributed multipole analysis of the isolated molecule charge density. Changes in structure were monitored via simulated X-ray diffraction patterns. The methodology was used to minimize the free energy at ambient conditions of a set of experimental and hypothetical 5-fluorouracil crystal structures, generated in a search for lattice-energy minima with the same model potential. Our results demonstrate that the majority (∼75%) of lattice-energy minima are thermally stable at ambient conditions, and hence, the free-energy (like the lattice-energy) surface is complex and highly undulating. Metadynamics trajectories (Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562) started from the free-energy minima only produced transitions that preserved the hydrogen-bonding motif, and thus, further developments are needed for this method to efficiently explore such free-energy surfaces. The existence of so many free-energy minima, with large barriers for the alteration of the hydrogen-bonding motif, is consistent with the range of motifs observed in crystal structures of 5-fluorouracil and other 5-substituted uracils.

1. Introduction The computational prediction of the organic crystal structure from the molecular diagram would have considerable impact on the development of fine chemicals.1 Current approaches are based on the assumption that the experimentally determined polymorphs correspond to minima at or near the global minimum on the lattice-energy surface. This approach, in conjunction with realistic models for the intermolecular interactions, has shown promise for certain categories of molecules2 but has also revealed that the number of distinct, energetically plausible lattice-energy minima typically exceeds the number of known polymorphs. In the majority of studies, the ranking of hypothetical crystal structures is solely based on the static lattice energy. The temperature dependence of the free energy of computationally generated polymorphs can be computed within the (quasi-)harmonic approximation.3-5 However, this method only probes the curvature of the static energy surface in the region around the minimum and ignores the real shape of the free-energy landscape, such as the barrier separating adjacent basins and their dependence on the thermodynamic conditions. A structure that lies in a shallow basin compared with kBT will be thermally unstable6 and under * To whom correspondence should be addressed. Tel: +44 (0)20 7679 4622. Fax: +44 (0)20 7679 7463. E-mail: [email protected]. † University College London. ‡ ETH. £ STFC Daresbury Laboratory.

thermal motion will spontaneously relax to an adjacent basin that may correspond to a fully structured crystal or a disordered state. Conventional molecular dynamics techniques can model spontaneous transitions when the free-energy barrier is small but provide no information on the structure of adjacent basins in phase space if the barrier is large. Hence, further transitions can only be observed by altering the thermodynamic conditions, for example, by heating or over-pressurization and thus causing an undesired change in the free-energy surface. The recently introduced method of metadynamics7 directly explores the free-energy surface at specified temperature and pressure within a genuinely dynamic simulation framework without being limited by the height of the free-energy barrier. Metadynamics comprises two simple building blocks, the identification of a set of reaction coordinates and the use of a history-dependent penalty potential on these coordinates. By driving the system out of even deep free-energy basins, it is possible to explore and reconstruct the free-energy landscape at the desired thermodynamic conditions. Metadynamics has been used in the modeling of processes involving the crossing of significant energy barriers (compared to kBT), such as reactions8,9 and solid-state phase transitions,10 which would otherwise require impractically lengthy simulation times to observe a spontaneous transformation to a nearby free-energy minimum. In the case of benzene, metadynamics has been shown to locate all experimentally identified polymorphs,11

10.1021/jp709764e CCC: $40.75 © 2008 American Chemical Society Published on Web 03/15/2008

Lattice-Energy Minima of 5-Fluorouracil

Figure 1. Hydrogen-bonding motifs encountered in the experimental and theoretically generated polymorphs of 5-fluorouracil. (a) Motif observed in form I. (b) Ribbon 1 observed in form II. (c) Ribbon 2. (d) Planar six-membered ring.

while no other hypothetical structures were generated despite the significantly larger number of lattice-energy minima.12 This paper expands the applicability of metadynamics in the case of the polymorphic system 5-fluorouracil as a representative example of a more typical organic molecule, whose solid-state exhibits hydrogen bonds and hence is not predominantly dispersion-bonded. Due to the strong electrostatic character of the intermolecular interactions, all molecular dynamics simulations were performed with an anisotropic intermolecular potential13 based on distributed multipoles computed from the isolated-molecule charge density.14 The accurate reproduction of the quantum mechanical molecular charge density with distributed multipoles was deemed necessary, as isotropic potentials based on atomic charges were found to underestimate the stability of hydrogen-bonded aggregates. By starting from static lattice-energy minima, we investigated the response of theoretically predicted and known polymorphs to thermal stress by minimizing their free energy and computing the curvature of the free-energy surface at the minimum as well as the thermal expansion coefficients. Metadynamics was then used to probe the structural relationship of adjacent basins in phase space by driving the system over the lowest free-energy saddle point. The onset of the induced phase transitions was identified by monitoring the X-ray powder pattern of the dynamically evolving system to detect abrupt structural changes. 2. Solid State of 5-Fluorouracil 5-Fluorouracil is known to exist in two crystalline forms. Form I15 crystallizes in space group P1h and exhibits an intricate, hydrogen-bonded sheet motif with four crystallographically independent molecules with regions of close F‚‚‚F contacts (Figure 1a). The recently discovered form II16 (P21/c, Z′ ) 1) adopts a hydrogen-bonded ribbon motif, with each molecule forming four NH‚‚‚O bonds with two adjacent molecules (hereafter ribbon 1, Figure 1b). A crystal structure prediction search for lattice-energy minima16 found form II as the most stable packing arrangement. The hydrogen-bonded ribbon of form II was also present in a large proportion of the hypothetical crystal structures within 10 kJ mol-1 of the global minimum. The next most dominant motif also has each molecule forming two hydrogen bonds to each of two neighboring molecules but differs in that the dimers share a common carbonyl acceptor

J. Phys. Chem. B, Vol. 112, No. 14, 2008 4299 (hereafter ribbon 2, Figure 1c), leaving one carbonyl group not hydrogen bonded. The ribbon 2 motif was found in the experimental structures of other uracils, including 5-chlorouracil,17 5-bromouracil,17 5-hydroxyuracil,18 thymine,19 and the 5-fluorouracil, 2,2,2-trifluoroethanol solvate.20 A third frequently encountered motif contains six-molecule rings arranged in planes, similar to the six-molecule ring found in the experimental crystal of uracil (63 sheet, Figure 1d).21 The diversity of hydrogen-bonded motifs in 5-substituted uracils increases the confidence in the predicted lattice-energy minima of 5-fluorouracil as thermodynamically plausible packing arrangements. Recent molecular dynamics simulations22 of 5-fluorouracil in water and dry and wet nitromethane suggest that the presence of water enhances the growth of form I, as both carbonyl and amine functional groups remain strongly solvated, and hence, their association into dimers is hindered. This rationalizes the experimental observation that form II was only crystallized from dry nitromethane.16 Thermal measurements16 showed that neither polymorph of 5-fluorouracil undergoes a phase transition prior to melting and that form I is more stable than form II. 3. Theoretical Methods 3.1. Free-Energy Minimization. All molecular dynamics runs were performed with supercells of approximately cubic shape comprising 32-48 molecules built to replicate the 60 most stable distinct lattice-energy-minimized structures16 so that all molecular dynamics supercell lengths were at least 15 Å. For consistency, both dynamical and static calculations used the same molecular model and intermolecular potential. In particular, the electrostatic contributions were evaluated from a set of atomic multipole moments up to hexadecapole computed via a distributed multipole analysis14 of the MP2/6-31G(d,p) molecular charge density. The repulsion-dispersion interactions were modeled with an empirical exp-6 potential (parameters given in Supporting Information) that was parametrized in conjunction with atomic multipole moments23 and supplemented with fluorine-fluorine interactions from an earlier parametrization.24 All molecular dynamics runs were performed by integrating the rigid-body equations of motion with 5-fluorouracil frozen in its MP2/6-31G(d,p)-optimized conformation with the package DL_MULTI,13,25,26 which can handle anisotropic intermolecular potentials. Electrostatic contributions were calculated with Ewald summation27 for all multipole ranks, while the van der Waals contributions were computed up to a 15 Å cutoff radius. DL_MULTI has been developed for the simulation of larger systems, for which the distance between opposite cell faces is at least twice the cutoff distance, which allows the use of periodic image convention28 to compute the interactions between atoms in neighboring cells. However, in the interest of computational efficiency and to reduce the formation of defects during metadynamics, we used relatively small cells. In order to allow the use of a reasonable cutoff distance that is of comparable size to cell lengths and to account for the appearance of obtuse or acute lattice angles during transitions, we modified DL_MULTI to dynamically alter the number of periodic replicas used in direct space summations, instead of relying on periodic image convention. The use of rigid bodies allowed good energy conservation with a 3 fs time step. The temperature was controlled with a Berendsen thermostat set at 0.1 ps. The first step of the investigation was to minimize the Gibbs free energy G(h), at 1 bar and 310 K, of the lattice-energy-

4300 J. Phys. Chem. B, Vol. 112, No. 14, 2008

Karamertzanis et al.

minimized structures. It can be shown10,28 that the gradient of the free energy with respect to the cell matrix h at temperature T and pressure P can be computed in a short molecular dynamics run in the canonical ensemble (NVT) as

∂G ) V[h-1(p - PI)]T ∂h

-

(1)

where p is the time-averaged pressure tensor. The importance of eq 1 is that it allows free-energy minimization by performing a series of NVT simulations with the cell changed with a steepest descent algorithm, using six degrees of freedom for a rotated cell such as it becomes upper triangular h h ≡ (h11,h22,h33,h12, h13,h23)T. The metadynamics cell changes were implemented by scaling the fractional atomic coordinates and placing the abinitio-optimized rigid molecules to maximize overlap with the distorted molecules produced by this scaling. The cell was equilibrated for 5 ps, followed by a 5 ps simulation time to compute the time-averaged pressure tensor. The starting steepest descent step was 0.1 Å, and it was halved each time two subsequent gradient vectors had an angle greater than 135°. Steepest descent relaxations required 10-40 steps to converge to an accuracy of 0.01 Å, that is, approximately 0.1% of the average cell dimension of the supercells considered. The steepest descent free-energy minimization performs gradual and controlled changes in the cell geometry, and hence, it is possible to converge to a shallow minimum that would be missed if greater cell variations were allowed. Moreover, the steepest descent can converge on the basis of negligible gradients when the structure lies in a lightly undulated, broad free-energy basin without providing any direct evidence on the shape of the basin. Hence, all free-energy-minimized structures were submitted to variable-cell isobaric isothermal molecular dynamics runs (NσT) for 100 ps under the same (ambient) thermodynamic conditions with a Berendsen barostat (0.3 ps). Only a small number of structures transformed to nearby basins, became defective, or exhibited pronounced oscillations in the cell geometry and were classified as thermally unstable. The thermal expansion coefficients of all thermally stable structures were used to probe their behavior under thermal stress and relate it to the computed lattice energy and packing arrangement. In principle, it is possible to integrate eq 1 along a trajectory connecting two basins to provide their free-energy difference. Unfortunately, our simulations showed that the thermodynamic integration is of comparable accuracy with the small stability difference of known and hypothetical polymorphs, that is, on the order of a few kJ mol-1. The accuracy of the integration could have been enhanced by reducing the step in cell changes and increasing the averaging time for the calculation of the pressure tensor, both increasing the computational expense significantly. This was beyond the aims of this study, and hence, although our minimization scheme can readily relax crystal structures, it has not been used to calculate relative free energies. 3.1.1. Thermal Expansion Coefficients. The minimization of the free energy with respect to the cell geometry was exploited to study the thermal expansion of the experimental and hypothetical crystal structures that remain stable at 310 K and 1 bar. The empirical repulsion-dispersion potentials used in this work were fitted to heats of sublimation and structural data at room or low temperatures and hence already contain some form of average thermal effects, which is reflected in the simulated cells. For this reason, we report the thermal expansion coefficients Rij that relate the resulting strain xij ) Rij∆T to the

temperature change ∆T, as the latter is less affected by the absolute density of the crystal compared with the cell lengths and angles. The thermal coefficients are reported in the IEEE standard setting,29 that is, for triclinic and monoclinic crystals z || c, y ⊥ ac plane [010] (for monoclinic y || unique axis b), and x lies on the [010] plane and forms a right-handed system. Axes are rearranged for orthorhombic systems so that z || c < x || a < y || b. The anisotropy of thermal expansion is correlated with the directionality of hydrogen bonds in the cell by visual inspection30 and, in the case of the experimentally determined 5-fluorouracil polymorphs, is studied as a function of temperature. The thermal expansion along an arbitrary direction was also linked to the temperature dependence of the simulated powder pattern by differentiating Bragg’s law, nλ ) 2d sin θ with respect to temperature to give

1 dd dθ ) -tan(θ) ) -tan(θ)Rθ dT d dT where Rθ is the thermal expansion in the direction normal to the plane that corresponds to diffraction angle θ. Hence, thermal expansion leads to shifting of the peaks to lower angles, with the shifting becoming increasingly pronounced for larger diffraction angles. 3.2. Free-Energy Hessian. For all thermally relaxed structures, we computed and diagonalized the Hessian matrix Hij ) ∂2G/∂hhi∂hhj to examine the curvature of the free-energy surface with respect to the cell geometry. Assuming that the free-energy basin is approximately quadratic close to the minimum, the freeenergy gradients can be approximated as linear functions of the = displacement from the equilibrium geometry (we use h to denote the cell geometry at the free-energy minimum). Due to the statistical nature of the pressure tensor, the second gradients ∂2G/∂hhi∂hhj cannot be readily estimated by numerical perturbations of the first gradients computed with eq 1. Our investigation showed that a more accurate procedure is to estimate the first derivatives ∂G/∂hhj at n points (where n is odd) =

hhkl ) hl + δil

2k δh n - 1 max

l ) 1, ..., 6

k ) -(n - 1)/2, ..., (n - 1)/2 and calculate the second gradients ∂2G/∂hhi∂hhj as the slopes of the least-squares fitted lines on these points. This procedure calculates all 36 elements of the Hessian matrix, which is then made symmetric, Hij ) 1/2(∂2G/∂hhi∂hhj + ∂2G/∂hhj∂hhi), to further reduce statistical noise. Preliminary investigations showed that the Hessian matrix is well converged (see Figure S1 in Supporting Information), with a maximum displacement δhmax ) 0.1 Å, n ) 5 and 45 ps average time for the pressure tensor. We test the hypothesis that structures exhibiting small Hessian eigenvalues can readily transform, which assumes that the curvature of the free-energy surface is linked to the freeenergy barrier for the transformation into the nearest basin. 3.3. Metadynamics. The final analysis of the free-energy surface was to investigate the structures that lie in the vicinity of the computed free-energy minima by inducing phase transitions into adjacent free-energy basins with metadynamics.10,11,31 All structures for which the smallest eigenvalue of the Hessian matrix was greater than a specified threshold (1.8 kJ mol-1 Å-2) were used as starting points for the metadynamics simulations. The rest of the minima (12 out of 60 distinct, static latticeenergy minima considered in this work) were rejected as metadynamics32 require a nonsingular Hessian matrix, that is, the curvature of the free-energy surface should not be negligible.


J. Phys. Chem. B, Vol. 112, No. 14, 2008 4301

In order to account for the significantly anisotropic curvature of the free surface, the system was evolved in the space of scaled = cell variables,32 s ) ΛD(h h - h), where Λ is a diagonal matrix containing the square root of the eigenvalues and D the orthonormal matrix with the eigenvectors of the Hessian at the = free-energy minimum h . New configurations were generated by a steepest-descent procedure sk+1 ) sk + δs(∂Gk/∂s/|∂Gk/ ∂s|), where

δGk/δs sk+1 ) sk + δs |δGk/δs|

(2)

is a modified free-energy surface. The second term augments the thermodynamic driving force with a history-dependent term that pushes the system out of the region of attraction of the minimum where it would otherwise be trapped and into a nearby one. The height W and width δs of the Gaussian functions of the history-dependent term were tailored to ensure that transitions are observed in a computationally manageable number of metadynamics steps, pushing the system smoothly over the lowest saddle point without abrupt cell changes.7 Preliminary investigations with a set of hypothetical 5-fluorouracil structures showed that δs ) 7.3 kJ1/2 mol-1/2 and W ) 6 kJ mol-1 ensure that the force exerted by the time-dependent potential never overwhelms that due to the free-energy landscape and therefore that the transition is driven by the free energy through the lowest energy barrier. Metadynamics simulations were carried on until the crystalline order was lost or for a maximum of 250 metasteps. Longer runs are not worthwhile because the scaling of the collective variables according to the free-energy curvature at the starting basin is no longer accurate due to significant structural modifications. The metadynamics transitions were examined to identify which structural variations were sampled. 3.3.1. Identification of Phase Transitions. Metadynamics simulations of silicon31 and silica32 showed that their phase transitions can readily be detected by monitoring the enthalpy and volume of the evolving system. This is due to the strong covalent character of all interactions in the lattice that leads to large differences in the density and packing energy of the various phases. In the case of organic compounds, intermolecular interactions are weak, and there are generally many distinct packing arrangements of similar density and energy. Hence, the enthalpy and volume do not always show appreciable changes in the course of metadynamics runs, and it is therefore not trivial to identify when a phase transition occurs. In this work, we detected transformations by comparing the time-averaged powder patterns (computed over a sufficiently large number of well-separated frames)33 of subsequent metadynamics steps. The similarity index between the powder patterns corresponding to the equilibrated structures in metadynamics steps k and k+1 was defined34 as

S

k,k+1

≡

∫ w(2θ)ck,k+1(2θ) d2θ ∫ w(2θ)ck,k(2θ) d2θ ∫ w(2θ)ck+1,k+1(2θ) d2θ

(3)

where cn,m(2θ) ≡ ∫ f n(2θ′)fm(2θ + 2θ′) d2θ′ is the powder pattern autocorrelation function. The weighting function w(2θ) ≡ max (1 - |2θ|/l,0) ensures that each point in the reference pattern k is compared with a neighborhood of the corresponding point in the target pattern k + 1. For all results reported in this paper, the width of this neighborhood, l, was set equal to 1° in 2θ, while the range of 2θ angles considered was from 10 to 50°. This similarity index varies from 0 (dissimilar) to 1

(identical) and does not need any prior normalization of the powder patterns. Once a transition is complete, the system oscillates around the equilibrium configuration of the new free-energy minimum, and hence, the similarity index between successive metadynamics steps is large. In contrast, immediately after surmounting a saddle point, the system swiftly relaxes toward the new freeenergy minimum, and hence, the similarity index drops. The passage through saddle points was also manifested by the alignment of the two components of the driving force, that is, the free energy and history-dependent contributions in eq 2. Metadynamics frames corresponding to new phases were quenched to 0 K with a quick NσT molecular dynamics run. The relaxed structures were compared against the static latticeenergy minima by computing the powder pattern and the rootmean-square difference between the optimally overlaid molecular coordination clusters.35 Every suspected new structure was thermally annealed for 50 ps at 100 K, followed by 50 ps at 300 K, and then cooled to 0 K in a series of NσT molecular dynamics simulations to overcome small energy barriers and reach the true equilibrium configuration. We opted for this annealing procedure and not a steepest-descent free-energy minimization because of the large number of relaxations involved and the qualitative character of this analysis. As a final step, the new structures were symmetry reduced with PLATON,36 and their thermodynamic stability contrasted with the hypothetical polymorphs generated in the search for static latticeenergy minima. The presence of a large number of molecules in the asymmetric unit indicated the formation of defects, such as stacking faults of parallel hydrogen-bonded sheets. 4. Results 4.1. Free-Energy Minimization of Experimentally Determined Polymorphs. Table 1 contrasts the predicted cell lengths and angles of both polymorphs of 5-fluorouracil with their experimental values as determined by single-crystal X-ray diffraction at 150 K and room temperature.15,16 Overall, the reproduction accuracy of the cell geometry is satisfactory during both static lattice-energy and dynamic free-energy minimizations, which suggests that our elaborate intermolecular potential models the strongly anisotropic interactions in the lattice sufficiently accurately. The overestimation in the cell volume at both temperatures can be attributed to the inclusion of average thermal effects in the repulsion-dispersion potential parametrization.23 In the case of form I, the error in density is approximately 11% and hence far exceeds the 3% density reduction due to thermal expansion between 150 K and room temperature. The thermal expansion of the experimentally determined 5-fluorouracil polymorphs I and II was studied by minimizing the free energy starting at 30 K and increasing to 310 K in 20 K steps (Figure 2 and Table 2). The predicted average volume expansivity β ) (1/VA)(VB - VA)/∆T for the temperature range of 150-310 K is 1.84 × 10-4 and 1.76 × 10-4 K-1 for forms I and II, respectively. These estimates agree reasonably well with the values of 2.33 × 10-4 and 1.64 × 10-4 K-1 derived from unit cell determinations at 150 K and room temperature.15,16 The larger expansivity of form I suggests that its melting point should be lower than that of form II, which is consistent with both its lower density and the predicted relative stability. However, previously reported16 thermal measurements showed the opposite monotropic relationship, as form I has both a higher melting point (onset of melting, form I, 282.6 °C; form II, 277.5 °C) and greater enthalpy of fusion (form I, 0.296 kJ g-1; form II, 0.200 kJ g-1).



Figure 2. Predicted thermal expansion coefficients as a function of temperature for (a) form I and (b) form II of 5-fluorouracil. Noise was removed by fitting the cell lengths and cell angles at the free-energy minima to quadratic polynomials of temperature that were subsequently used to compute the thermal expansion coefficients. Thermal expansion coefficients were computed from the fitted cell geometries between two successive temperature points (0, 30, 50, ..., 310 K) and reported for the midpoint temperature.

TABLE 1: Reproduction of Cell Geometry for Forms I and II of 5-Fluorouracil by Lattice-Energy (LE) and Free-Energy (FE) Minimization cell geometry a (Å)

b (Å)

c (Å)

0 K, 0 bar 150 K, 1 bar 150 K, 1 bar RT, 1 bar 310 K, 1 bar

predicted (LE) experimental16 predicted (FE) experimental16 predicted (FE)

8.839 8.633 8.939 8.786 9.066

Form I, P1h, Z′ ) 4 9.279 13.048 9.156 12.580 9.346 13.067 9.220 12.660 9.440 13.112

0 K, 0 bar 150 K, 1 bar 150 K, 1 bar RT, 1 bar 310 K, 1 bar

predicted (LE) experimental16 predicted (FE) experimental15 predicted (FE)

5.350 5.043 5.427 5.154 5.556

Form II, P21c, Z′ ) 1 15.262 6.508 14.935 6.605 15.307 6.571 15.001 6.654 15.388 6.640

a

R (°)

β (°)

γ (°)

Va (Å3)

83.39 80.88 83.59 81.40 83.79

82.86 79.98 83.16 80.53 83.44

91.69 90.02 91.08 89.41 90.15

131.73 120.80 134.57 125.01 138.53

110.29 108.88 111.21 110.36 112.83

124.60 117.67 127.22 120.58 130.81

Cell volume per molecule.

TABLE 2: Reproduction of Thermal Expansion Coefficients for Forms I and II of 5-Fluorouracil average thermal expansion coefficients (10-4 K-1)

c

R11

R22

R33

experimentalb predictedc

1.29 0.93

0.58 0.68

experimentalb predictedc

0.84 0.76

0.29 0.33

R12

R13

R23

βa (10-4 K-1)

Form I, P1h, Z′ ) 415,16 0.43 0.48 0.22 0.55

-0.26 -0.11

-0.31 -0.09

2.33 1.84

Form II, P21c, Z′ ) 116 0.49 0.66

-1.04 -1.06

1.65 1.76

a Thermal expansivity. b Computed over the temperature range of 150-300 K, taking room temperature as 300 K for the experimental investigations. Computed over the temperature range of 150-310 K.

The data in Table 2 show that the agreement is reasonable for all thermal expansion coefficients Rij that do not vanish due to symmetry, despite the significant errors in the predicted density. The larger errors for form I can be attributed to the greater influence of the close fluorine-fluorine interactions that were not sufficiently sampled in the force-field parametrization.24 Thermal expansion is very anisotropic, with some components of the second-order tensor being more temperature dependent than others. The symmetry-allowed shear strains (offdiagonal thermal expansion coefficients) appear comparable in magnitude to the positive tensile expansions, suggesting a strong dependence of the cell shape on temperature. The predicted thermal expansion of form I along crystallographic axis c is only R33 ≈ 0.22 × 10-4 K-1, which is expected as cell deformations in this direction require changes in the hydrogen-

bonded sheet motif which lies on plane (1 1 0). Heating causes significantly stronger expansion in the direction normal to this plane (1.33 × 10-4 K-1). It is also consistent that thermal expansion of form II in the direction normal to the best-fit plane (12 2 -24) passing through the sheets formed by the hydrogenbonded ribbons is 1.14 × 10-4 K-1 and hence much larger than thermal expansion in the direction of the ribbons, which is only 0.10 × 10-4 K-1. The effect of thermal expansion on the time-averaged simulated powder pattern for forms I and II of 5-fluorouracil is shown in Figure 3. The strong anisotropic character of thermal expansion results in neighboring peaks that shift at different rates with temperature. The resolution of peaks that overlap at a different temperature can be readily predicted with the proposed free-energy minimization scheme to the accuracy



Figure 3. Effect of temperature on the simulated powder pattern of 5-fluorouracil (a) form I and (b) form II. The inset shows two Miller planes for form II that are associated with markedly different thermal expansion in the direction orthogonal to the planes, which explains the variation in the temperature shifts of the corresponding powder X-ray diffraction peaks.

TABLE 3: Thermodynamic Stability, Thermal Expansivity, and Curvature of the Free-Energy Surface with Respect to the Cell Matrix for the Experimental and 20 Most Stable Hypothetical Crystal Structures of 5-Fluorouracil 0 K, 0 bar structurea

space group

form I

P1h, Z′ ) 4

form II

P21/c

1. am56 2. ca13 3. am108 4. am93 5. ca33 6. ak24 7. ak90 8. am51 9. aq61 10. fa47 11. am36

P21/c P1h P21/c P21/c P1h P21/c P21/c P21/c P212121 P21/c P21/c

12. dc71 13. aq90 14. fa20 15. fc19 16. am64

C2/c P212121 P21/c P21/c P21/c

17. ab94f 18. ab98 19. ab13 20. ak39

P1h P1h P1h P21/c

packing motif planar sheet, F-F contacts ribbon 1 ribbon 1 ribbon 2 ribbon 1 ribbon 1 ribbon 2 ribbon 2 ribbon 2 ribbon 1 ribbon 2 ribbon 2 planar six-member ring ribbon 2 ribbon 2 distorted ribbon 2 ribbon 2 planar six-member ring distorted ribbon 2 distorted ribbon 2 ribbon 2 ribbon 2

310 K, 1 bar

Vb (Å3)

lattice energy (kJ mol-1)

Vb (Å3)

lattice enthalpyc (kJ mol-1)

βd (10-4 K-1)

131.7

Experimental -96.58

138.6

-81.90

1.68

3.68

124.6

-102.48

130.8

-86.80

1.60

3.85

123.8 123.1 125.8 126.6 126.8 126.6 127.3 124.4 127.5 128.7 124.4

Hypothetical -102.40 -100.50 -100.34 -100.23 -100.19 -99.64 -99.26 -99.11 -98.99 -98.96 -98.78

-87.00 transformed to ca33 -84.65 -84.87 -85.02 thermally unstable thermally unstable -83.81 -83.48 thermally unstable -83.48

1.60

2.06 3.67 12.09 4.31 0.36 0.88 1.97 21.84 11.30 3.14 7.45

128.5 128.4 125.7 128.2 122.7

-98.75 -98.60 -98.53 -98.49 -98.35

1.65

132.7 134.3 129.3

-82.53 thermally unstable -82.56 -83.03 -83.02

125.5 125.8 125.6 126.8

-98.30 -97.98 -97.89 -97.83

132.7 132.1 132.9

transformed to ab98 -82.49 -82.44 -82.50

129.9 132.1 132.7 132.6 131.1 134.3 130.9 135.1

1.60 1.55 1.46 1.72 1.72 1.67

min λe (kJ mol-1 Å-2)

1.81 1.52 1.74

16.60 -2.14 4.71 0.08 4.44

1.77 1.67 1.55

6.73 10.11 3.79 5.89

a Lattice-energy rank and structure name as in original CSP search.16 b Cell volume per molecule. c Average enthalpy at the thermally expanded cell over a 5 ps period. d Average volume expansivity over the temperature range of 0-310 K. e Smallest eigenvalue of the Gibbs free-energy Hessian with respect to the nonzero components of the cell matrix at the free-energy minimum. f Transition occurred after 60 ps. The proximity of the two basins on the free-energy surface is comparable to the perturbation in cell geometry to numerically compute the Hessian, which may explain the large λ.

allowed by the intermolecular potential. For example, the powder pattern of the lattice-energy-minimized form II (Figure 3) exhibits two characteristic peaks that correspond to Miller planes (1 1 1) and (1 0 -2), which are well resolved at 0 K but overlap partially when simulated at 310 K. The near invariance of the (1 1 1) peak with temperature reflects the insignificant expansion of the hydrogen bonds, whereas the significant shift of the (1 0 -2) reflection reflects the expansion approximately orthogonal to the hydrogen-bonded sheets.

4.2. Thermodynamic Stability of CSP Structures under Ambient Conditions. Table 3 shows the results of the freeenergy minimization at 310 K and 1 bar of the experimental and 20 most stable hypothetical crystal structures of 5-fluorouracil; Table S1 in Supporting Information completes this list for all 60 lattice-energy minima within 8 kJ mol-1 of the global minimum. All but two steepest-descent minimizations of the free energy terminated successfully on the basis of negligible gradients. However, a significant proportion of the resulting


Figure 4. Cell length deformation (defined as the root-mean-square relative change in the cell lengths) during the NσT simulations starting from the free-energy minima as a function of the smallest eigenvalue of the free-energy Hessian. All structures that showed cell length deformation greater than 2.5% or a standard deviation of any cell length or angle in excess of 0.8 Å and 5°, respectively, during the NσT run were considered thermally unstable. Filled triangles denote thermally unstable structures for which we could identify the resulting structure in the list of hypothetical polymorphs.

structures (∼25%) corresponded to shallow free-energy basins and led to significant geometry relaxation (root-mean-square deformation of lattice lengths greater than 2.5%) when subjected to NσT molecular dynamics at 310 K and 1 bar due to transformation to other crystalline phases or amorphization. In some instances, the cell lengths and angles showed large oscillations, suggesting the presence of a broad free-energy basin; these structures were also classified as thermally unstable if the standard deviation for any of the lattice lengths of the approximately cubic 15 Å starting cell was greater than 0.8 Å or the standard deviation of a lattice angle exceeded 5°. While, for the majority of structures, the simulated powder pattern during the NσT run was the same as that of the steepest-descent free-energy minimum, for the thermally unstable structures, we observed new peaks and significant peak broadening (see Figure S2 Supporting Information). This is consistent with the steepestdescent minimization having arrived in a very shallow and probably slightly undulating basin. Although most theoretical polymorphs remain stable under thermal motion at ambient conditions, the curvature of the freeenergy surface at the corresponding free-energy minima varies significantly. The smallest Hessian eigenvalue gives the curvature of the free-energy surface along the direction in which cell deformation is most likely to occur and hence provides insights into the susceptibility of the structure to deformation. The largest eigenvalue, which was found to correspond to cell geometry changes associated with significant disruption of the hydrogen-bonding pattern, is typically two orders of magnitude greater than the smallest that is associated with energetically less demanding packing rearrangements. Most of the thermally unstable structures within 8 kJ mol-1 of the global lattice-energy minimum exhibit a ribbon 2 hydrogen-bond motif (see Table 3 and S1 in Supporting Information), with the smallest Hessian eigenvalue corresponding to a facile sliding of hydrogen-bonded sheets. There is a correlation between the smallest eigenvalue of the Hessian and the root-mean-square cell length deformation during the NσT runs (Figure 4). In particular, the structures with a large eigenvalue (greater than 8 kJ mol-1 Å-2) are always stable during the NσT runs, while those with a very small eigenvalue often lead to significant structural relaxations during the NσT simulation. However, there is an intermediate range of eigen-

Karamertzanis et al. values which corresponds to both thermally stable and unstable structures, which includes both experimental polymorphs. For a small number of structures, the Hessian matrix is not positive definite despite the successful termination of the steepest descent minimization. This serves as warning evidence that, despite the elaborate numerical evaluation of the second freeenergy gradients, the accuracy of the Hessian matrix remains limited, especially for low curvature minima that belong in freeenergy surface regions that contain shallow basins in close proximity. A closer inspection showed that in such cases, the assumed linear dependence of the first derivatives on the perturbation step is not accurate. Unfortunately, reducing the perturbation step would increase the statistical noise and would require an impractically long simulation time. We needed to confirm that the transformation of thermally unstable structures into more stable configurations was not the onset of melting due to inadequacies in our simulation model. A qualitative investigation by heating the two experimental structures in a series of free-energy minimizations with 40 K temperature increments showed no evidence of thermal disorder below 550 K. We did not attempt to heat the structures a further 300 K in order to compare with the experimental melting points,16 given the technical challenges in modeling solid-liquid coexistence for molecular systems.37 However, the thermal stability of the known structures up to 550 K strongly suggests that the observed phase transitions or loss of symmetry at 310 K for several of the hypothetical polymorphs is not due to melting but shows thermally instability. 4.3. Phase Transitions with Metadynamics. Metadynamics runs starting from the experimental structures did not show any transition to other phases. Form I appears isolated in phase space, and after a few tens of metasteps, it showed extended defects and eventually became disordered. Form II was produced during a metadynamics run starting from the second most stable predicted structure (am56; see Table 3), but the backward transition was never observed, suggesting that form II belongs to a deep minimum. This is also confirmed by the fact that independent metadynamics runs with different simulation cells only showed transitions between different supercells that all corresponded to form II. The metadynamics trajectory starting from the hypothetical crystal structure ab94 (Figure 5) is characteristic of the structural transformations observed between 5-fluorouracil crystals. The starting structure contains deformed, nonplanar type-2 ribbons (Figure 1). Adjacent ribbons are connected with a strained hydrogen bond between a bifurcated amine hydrogen and the carbonyl acceptor that is typically unused in the planar ribbon 2 motif. The significant changes in lattice lengths and angles during metadynamics are linked to commensurate packing rearrangements, which, however, barely vary the hydrogenbonded ribbon motif. The first transition into the structurally very similar ab98 theoretical polymorph does not cause any change in the hydrogen-bond motif. The transition to the 2 kJ mol-1 more stable ca33 structure breaks the strained hydrogen bond to the bifurcated amine hydrogen. The energy penalty for the breaking of this hydrogen bond is compensated by the adoption of a planar, and assumingly energetically favorable, geometry by the hydrogen-bonded ribbons. The continuation of the metadynamics trajectory, despite a significant sliding and modulation of the sheets on which the ribbons lie, leads to no appreciable change in the ribbons and to little change in the overall packing. The structures connected by metadynamics trajectories are shown in Figure 6. Structures with hydrogen bonding in two or



Figure 6. Pictorial representation of the metadynamics results, with the hypothetical and observed structures grouped by hydrogen-bond motif. Structures that interconvert to each other during metadynamics runs are contained in rectangular boxes. New crystalline structures found after quenching selected metadynamics steps are shown in bold (the name gives the starting structure, and the number after the underscore gives the metadynamics step that led to it). Structures that were not crystalline or that had many stacking faults are not shown.

divided into regions containing structures exhibiting similar hydrogen-bonding patterns. 5. Discussion Figure 5. Metadynamics trajectory of structure ab94, showing changes in the X-ray powder pattern (top panel) between successive metadynamics structures. The continuous line serves to guide the eye to the trend in the X-ray similarity coefficient between subsequent metadynamics steps. The vertical arrows indicate the metadynamics structures that were quenched to 0 K and the resulting structures. The multiple transitions to the same structure are due to transformations into crystallographically equivalent cells. The bottom four panels show the difference in lattice energy, unit cell parameters, and volume per molecule of the structure at each metadynamics step relative to the values for the free-energy minimized ab94 structure. When metadynamics causes a transition between structural basins, there are large changes in unit cell lengths and angles (not necessarily associated with equally significant variations in the lattice energy and volume per molecule), whereas subsequent metadynamics steps within the same basin lead to gradual structural variations.

three dimensions, for example, the six-membered ring in Figure 1d or a 3D hydrogen-bonded network, are resistant to transformations, and eventually, when the time-dependent potential overwhelms the packing forces, defects develop, and the crystalline order is abruptly lost. On the other hand, structures containing hydrogen-bonded ribbons arranged in parallel sheets transform more readily via a sliding of the sheets, especially for ribbon 2 structures. For these structures, a few metadynamics steps were enough to observe a transition. Quenching of selected metadynamics frames (that correspond to similarity index minima) to 0 K and analyzing the symmetry of the resulting structures revealed that most new structures contain a large number of molecules in the asymmetric unit due to stacking faults that do not lead to a significant increase in lattice energy. This behavior is consistent with the analysis of the Hessian eigenvectors, which always indicates the sliding of hydrogenbond sheets as the energetically most plausible deformation mechanism. Overall, we observed no clear transition between distinct hydrogen-bonded motifs, such as transformation of ribbon 1 to ribbon 2. Hence, the phase space appears to be

5.1. Free-Energy Minimization. The dependence of the Gibbs free-energy gradient with respect to cell geometry on the readily calculated average pressure tensor is exploited in developing a reliable steepest-descent, free-energy minimization algorithm. The algorithm is able to model thermal expansion in the entire temperature range in which a crystal structure is stable, that is, is valid beyond the low-temperature region of harmonic molecular librations. The encouraging agreement of the predicted thermal expansion coefficients for both 5-fluorouracil polymorphs shows that the anisotropy in thermal expansion due to strong, directional intermolecular interactions is accurately modeled. One practically important application of this methodology may be to assist structure determination from powder X-ray data, when peak overlap prevents indexing and it is not possible to obtain powder patterns at a sufficient range of temperatures to detect the peak overlap experimentally. When indexing is not possible, Rietveld refinement38 can be started from the theoretically predicted lattice-energy minimum that closest matches the measured X-ray powder pattern.39 Simulating the thermal variation in the powder pattern could improve the confidence in matching the experimental with the corresponding simulated powder X-ray pattern and additionally facilitate Rietveld refinement, if the latter is started from the thermally expanded cell at the measurement temperature. 5.2. Thermal Stability of Hypothetical Polymorphs. The application of metadynamics in the case of benzene11 showed that the number of free-energy basins is considerably smaller than the number of lattice-energy minima.12 Hence, it was argued6 that the complexity of the lattice-energy surface2 is due to the search-locating minima that lie in such shallow basins, separated from more stable configurations by small barriers readily surmounted by thermal motion. Our results demonstrate that almost all 5-fluorouracil lattice-energy minima within 8 kJ mol-1 of the global minimum correspond to a unique free-energy basin. However, approximately 25% of these are shallow and

4306 J. Phys. Chem. B, Vol. 112, No. 14, 2008 transform to more stable configurations during a short variablecell isobaric-isothermal molecular dynamics simulation. Ribbon 2 structures are more likely to be thermally unstable due to the facile sliding of the hydrogen-bonded sheets. On the other hand, there are no structures that contain the experimentally observed ribbon 1 motif that either transform or lose their crystalline order. Thus, in marked contrast to benzene, the consideration of thermal effects does not appreciably simplify the energy landscape of crystalline 5-fluorouracil. A wide range of packings of different hydrogen-bonded motifs are free-energy minima, and therefore, another explanation for why more polymorphs have not been observed is required. The curvature of the Gibbs free-energy surface can serve as useful a posteriori analysis of crystal structure prediction studies to eliminate a proportion of the almost equienergetic latticeenergy minima as susceptible to phase transformations under thermal motion. Those 5-fluorouracil structures for which all eigenvalues of the free-energy Hessian are large (greater than 8 kJ mol-1 Å-2) do not relax under thermal motion (Figure 4). Hence, a steep basin at the vicinity of the minimum is associated with a barrier that is sufficiently large to prevent thermal instability. A large proportion of free-energy minima that belong to wells of small curvature correspond to thermally unstable structures in the NσT simulation, implying that they are shallow minima in a slightly undulating free-energy surface region. However, there is a region of intermediate curvatures for which thermal stability cannot be predicted by diagonalizing the freeenergy Hessian (Figure 4), which includes the experimental structures. The thermal analysis for 5-fluorouracil structures appears to be complex because the free-energy surface contains slightly undulating regions with shallow minima as well as a large number of deep basins. It has been suggested that simulating the melting of crystal structures could be used to distinguish real from hypothetical polymorphs, as the melting profiles of four hypothetical structures of coumarin40 were observed to start at significantly lower temperatures than the experimental structure, despite similar lattice energies.41 For the computational polymorphs that remain thermally stable at 310 K and 1 bar, we could assume that the average thermal expansion in the temperature range of 0-310 K is inversely proportional to the melting point. The average thermal expansivity (Tables 3 and S1 in Supporting Information) of the two experimentally determined polymorphs (form I, 1.68 × 10-4 K-1; form II, 1.60 × 10-4 K-1) lies in the lower part of the 1.4-2.4 × 10-4 K-1 range that contains the majority of the theoretical polymorphs within 8 kJ mol-1 from the global lattice-energy minimum. However, it seems improbable that melting simulations would clearly distinguish the observed polymorphs of 5-fluorouracil from the hypothetical structures. Overall, the data on the simulated polymorphs of 5-fluorouracil (Tables 3 and S1) show that many crystal energyproperty correlations that apply to dispersion-bound systems, like benzene11 or coumarin,40 are less applicable to 5-fluorouracil. 5.3. Induced Phase Transitions with Metadynamics. The hydrogen-bonding pattern remains largely invariant during the phase transitions and structural rearrangements of thermally unstable structures or during the metadynamics trajectories of thermally stable computational polymorphs. There is no unambiguous breaking and forming of hydrogen bonds, with, at most, subtle modifications of pre-existing but highly deformed ribbons into the more favorable planar geometry. This is in marked contrast to the successful exploration of the free-energy landscape of benzene,11 for which all stable free basins were


Figure 7. Energy penalty (kJ mol-1 of the three molecule cluster) to deform the hydrogen bonds of a ribbon 1 cluster of 5-fluorouracil; (a) the elongation of one of the two sets of parallel N-H‚‚‚O hydrogen bonds and (b) the tilting of the end molecule of the ribbon with respect the plane of the other two. The calculations are carried out using ORIENT,45 with the same model intermolecular potential as that in the MD simulations. To illustrate the difference between the distributed multipole (Rank 4) electrostatic model and an atomic charge (Rank 0) representation of the same MP2-6-31G(d,p) charge density, the deformation energies with the CHELPG46 atomic charges are also shown.

accessible in a single metadynamics trajectory. The inability of metadynamics to break and form hydrogen bonds (Figure 6) suggests a partitioning of the free-energy landscape into regions of similar hydrogen-bonding motifs. The resilience of the 5-fluorouracil hydrogen-bond geometries to structural relaxation is explained (Figure 7) by the intermolecular energy penalty for deformation of the ribbon 1 hydrogen-bonding motif within a three molecule cluster. Elongation or bending of one of the double hydrogen bonds by only 0.5 Å or 55°, respectively, decreases the stability of the three molecule cluster by 15 kJ mol-1. This energy penalty for the hydrogen-bond deformation is significant given that the lattice-energy range for all 60 structures presented in this work is only 8 kJ mol-1. (Inciden-

Lattice-Energy Minima of 5-Fluorouracil tally, the less realistic atomic charge model underestimates the deformation energy barrier and particularly the bending of the ribbons.) The thermal stability of the hydrogen-bonded motifs is consistent with both 5-fluorouracil polymorphs showing no transition in the whole temperature range up to their melting point. It is also consistent with the hydrogen-bonding motif being determined early in the nucleation process, as suggested by simulations22 of 5-fluorouracil in the solvents that produce the different polymorphs. The range of packing rearrangements seen in the metadynamics trajectories (Figure 5) with small energy changes and little variation in the hydrogen-bonding motif are likely to be the type that can most readily occur within large nucleation clusters. The very different result of molecular dynamics studies on the real and hypothetical crystal structures of 5-fluorouracil and benzene or coumarin can be explained in terms of the strength of the hydrogen bonding. It is important to note that this is specific to hydrogen bonds of the strength and directionality of 5-fluorouracil (Figure 7), and qualitatively different results are expected for less strongly associating systems. Facile transformations between hydrogen-bonding motifs have been simulated in molten acetic acid,42 albeit with an atomic charge electrostatic model. The ease of transformation between possible hydrogenbonding motifs due to unusually weak N-H‚‚‚OdC hydrogen bonds explains both the occurrence of a plastic phase for 3-azabicyclo(3.3.1)nonane-2,4-dione and the lack of a polymorph with the alternative hydrogen-bonding motif.43 5.4. Future Metadynamics Developments. The presence of two different energy scales for the intermolecular interactions in the organic solid state is a major limitation for the applicability of metadynamics to efficiently explore phase space. The freeenergy barrier separating structures that only differ in the packing of the same hydrogen-bond motif is considerably smaller compared with the barrier separating structures with different motifs. Hence, if the Gaussians’ height (W in eq 2) is suitable to reconstruct the small (dispersion-related) modulations of the free-energy surface, it will not be possible to observe changes in hydrogen bonding in a computationally feasible number of metadynamics steps. On the other hand, choosing tall Gaussians will be detrimental for the ability of the method to visit wells separated by moderate barriers, which may correspond to different polymorphs. Hence, further methodological developments are necessary. One possible route is to investigate the use of additional collective variables to reflect the hydrogen-bonding patterns in the lattice. Alternatively, metadynamics may be combined with replica exchange methods44 by running parallel simulations with progressively weakened hydrogen-bonding interactions in order to produce one replica in which all intermolecular interactions lie on the same energy scale, which will allow all free-energy basins to be visited. These future developments will enhance the ability of metadynamics to locate all stable free-energy basins at the thermodynamic conditions of interest but will not change their number. Strongly associating molecules, such as 5-fluorouracil, are likely to have a more complex free-energy surface compared with dispersion-bound systems, such as benzene,11 with the majority of lattice-energy minima belonging to sufficiently deep free-energy basins that they potentially correspond to long-lived, metastable polymorphs. 6. Conclusions In this paper, we reported the development of a free-energy minimization protocol to thermally relax crystal structures

J. Phys. Chem. B, Vol. 112, No. 14, 2008 4307 exhibiting strong, anisotropic intermolecular interactions due to hydrogen bonding. The methodology reproduces well the structure and thermal expansion coefficients of polymorphs exhibiting very different packing arrangements. When combined with crystal structure prediction methodologies, the modeling of thermal expansion can assist the determination of crystal structures from nonindexable X-ray powder diffraction data by establishing the variable temperature dependence of peaks due to the anisotropy of interactions in the lattice. The free-energy surface in the vicinity of known and hypothetical polymorphs of 5-fluorouracil was examined by studying the curvature of the free-energy surface with respect to cell changes and additionally by driving the system into nearby basins with metadynamics. Our results suggest that the free-energy surface of 5-fluorouracil is complex, with a large number of basins separated with large or small barriers, depending on the similarity of the corresponding hydrogen-bond motifs. The majority of static lattice-energy minima belong to distinct free-energy wells, although a quarter of these were sufficiently shallow to show thermal instability. Acknowledgment. Dr. Ashley Hulme and Dr. Maryjane Tremayne are thanked for useful discussions, and Dr. Louise Price is thanked for assistance in preparing the manuscript. Funding by the Basic Technology program of the Research Councils U.K. (CPOSS project www.cposs.org.uk) and the Swiss National Science Foundation (Grant Number 200021109673) is also gratefully acknowledged. This work made use of the facilities of CSCS (Swiss National Supercomputing Centre) and HPCx (U.K.’s national high-performance computing service provided by EPCC at the University of Edinburgh and by CCLRC Daresbury Laboratory). Supporting Information Available: Figure showing numerical tests to accurately compute the free-energy Hessian matrix. Simulated powder patterns for a subset of the thermally unstable theoretical and experimental polymorphs. A continuation of Table 3 of the manuscript showing all structures within 8 kJ mol-1 of the global lattice-energy minimum. Table with cell geometries for the static-lattice-energy-minimized (0 K), free-energy-minimized (310 K), and NσT-relaxed (310 K) crystal structures. Table with average thermal expansion coefficients for all hypothetical polymorphs in the temperature range of 0-310 K. Table with repulsion-dispersion potential parameters. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Day, G. M. et al. Acta Crystallogr., Sect B 2005, 61, 511. (2) Price, S. L. AdV. Drug DeliVery ReV. 2004, 56, 301. (3) Van Eijck, B. P. J. Comput. Chem. 2001, 22, 816. (4) Day, G. M.; Price, S. L.; Leslie, M. J. Phys. Chem. B 2003, 107, 10919. (5) Gavezzotti, A.; Filippini, G. J. Am. Chem. Soc. 1995, 117, 12299. (6) Mooij, W. T. M.; Van Eijck, B. P.; Price, S. L.; Verwer, P.; Kroon, J. J. Comput. Chem. 1998, 19, 459. (7) Laio, A.; Parrinello, M. P. Natl. Acad. Sci. U.S.A. 2002, 99, 12562. (8) Ensing, B.; Laio, A.; Gervasio, F. L.; Parrinello, M.; Klein, M. L. J. Am. Chem. Soc. 2004, 126, 9492. (9) Stirling, A.; Iannuzzi, M.; Laio, A.; Parrinello, M. ChemPhysChem 2004, 5, 1568. (10) Martonak, R.; Laio, A.; Bernasconi, M.; Ceriani, C.; Raiteri, P.; Zipoli, F.; Parrinello, M. Z. Kristallogr. 2005, 220, 489. (11) Raiteri, P.; Martonak, R.; Parrinello, M. Angew. Chem., Int. Ed. 2005, 44, 3769. (12) Van Eijck, B. P.; Spek, A. L.; Mooij, W. T. M.; Kroon, J. Acta Crystallogr., Sect. B 1998, B54, 291.

4308 J. Phys. Chem. B, Vol. 112, No. 14, 2008 (13) Gray, A. E.; Day, G. M.; Leslie, M.; Price, S. L. Mol. Phys. 2004, 102, 1067. (14) Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047. (15) Fallon, L. Acta Crystallogr., Sect. B 1973, B 29, 2549. (16) Hulme, A. T.; Price, S. L.; Tocher, D. A. J. Am. Chem. Soc. 2005, 127, 1116. (17) Sternglanz, H.; Bugg, C. E. Biochim. Biophys. Acta 1975, 378, 1. (18) Copley, R. C. B.; Deprez, L. S.; Lewis, T. C.; Price, S. L. Cryst. Eng. Commun. 2005, 7, 421. (19) Portalone, G.; Bencivenni, L.; Colapietro, M.; Pieretti, A.; Ramondo, F. Acta Chem. Scand. 1999, 53, 57. (20) Hulme, A. T.; Tocher, D. A. Acta Crystallogr., Sect. E 2005, 61, O3661. (21) Stewart, R.; Jensen, L. Acta Crystallogr. 1967, 23, 1102. (22) Hamad, S.; Moon, C.; Catlow, C. R. A.; Hulme, A. T.; Price, S. L. J. Phys. Chem. B 2006, 110, 3323. (23) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J. Phys. Chem. 1996, 100, 7352. (24) Williams, D. E.; Houpt, D. J. Acta Crystallogr., Sect. B 1986, 42, 286. (25) Leslie, M. Mol. Phys. 2008, submitted. (26) Liem, S. Y.; Popelier, P. L. A.; Leslie, M. Int. J. Quantum Chem. 2004, 99, 685. (27) Ewald, P. Ann. Phys. 1921, 64, 253. (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (29) Newnham, R. E. Properties of Materials: Anisotropy, Symmetry, Structure; Oxford University Press: New York, 2004. (30) Macrae, C. F.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Shields, G. P.; Taylor, R.; Towler, M.; De Streek, J. J. Appl. Crystallogr. 2006, 39, 453. (31) Martonak, R.; Laio, A.; Parrinello, M. Phys. ReV. Lett. 2003, 90, 075503. (32) Martonak, R.; Donadio, D.; Oganov, A. R.; Parrinello, M. Nat. Mater. 2006, 5, 623.

Karamertzanis et al. (33) Torrisi, A.; Leech, C. K.; Shankland, K.; David, W. I. F.; Ibberson, R. M.; Benet-Buchholz, J.; Boese, R.; Leslie, M.; Catlow, C. R. A.; Price, S. L. J. Phys. Chem. B 2008, DOI: 10.1021/jp710017y. (34) De Gelder, R.; Wehrens, R.; Hageman, J. A. J. Comput. Chem. 2001, 22, 273. (35) Chisholm, J. A.; Motherwell, S. J. Appl. Crystallogr. 2005, 38, 228. (36) Spek, A. L. PLATON, A Multipurpose Crystallographic Tool; Utrecht University: Utrecht, The Netherlands, 2003. (37) Eike, D. M.; Maginn, E. J. J. Chem. Phys. 2006, 124, 164503. (38) Rietveld, H. M. J. Appl. Crystallogr. 1969, 2, 65. (39) Harris, K. D. M.; Tremayne, M.; Kariuki, B. M. Angew. Chem., Int. Ed. 2001, 40, 1626. (40) Gavezzotti, A. J. Am. Chem. Soc. 2000, 122, 10724. (41) If the entropy of melting of all 5-fluorouracil computational polymorphs is assumed to be equal to the average value of the two experimental forms16 (form I, 69.3 J K-1 mol-1; form II, 47.0 J K-1 mol-1), the difference in their melting point should be on the order of 17 K for each 1 kJ mol-1 of enthalpy difference. This suggests that the 8 kJ mol-1 lattice-energy range considered in this work corresponds to a melting point range of approximately 140 K, which is expected to be detectable in molecular dynamics simulations, despite the approximations involved. (42) Gavezzotti, A. J. Mol. Struct. 1999, 486, 485. (43) Hulme, A. T.; Johnston, A.; Florence, A. J.; Fernandes, P.; Shankland, K.; Bedford, C. T.; Welch, G. W. A.; Sadiq, G.; Haynes, D. A.; Motherwell, W. D. S.; Tocher, D. A.; Price, S. L. J. Am. Chem. Soc. 2007, 129, 3649. (44) Earl, D. J.; Deem, M. W. Phys. Chem. Chem. Phys. 2005, 7, 3910. (45) Stone, A. J.; Dullweber, A.; Engkvist, O.; Fraschini, E.; Hodges, M. P.; Meredith, A. W.; Nutt, D. R.; Popelier, P. L. A.; Wales, D. J. ORIENT: A Program for Studying Interactions between Molecules, version 4.6; University of Cambridge: Cambridge, U.K., 2007. (46) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361.

The Thermal Stability of Lattice-Energy Minima of 5-Fluorouracil

Recommend Documents