Ctirad
Červinka
a and
Gregory J. O.
Beran
*b
aDepartment of Physical Chemistry, University of Chemistry and Technology Prague, Technická 5, CZ-16628 Prague 6, Czech Republic. E-mail: cervinkc@vscht.cz
bDepartment of Chemistry, University of California, Riverside, California 92521, USA. E-mail: gregory.beran@ucr.edu
First published on 16th April 2018
Organic crystals frequently adopt multiple distinct polymorphs exhibiting different properties. The ability to predict not only what crystal forms might occur, but under what experimental thermodynamic conditions those polymorphs are stable would be immensely valuable to the pharmaceutical industry and others. Starting only from knowledge of the experimental crystal structures, this study successfully predicts the methanol crystal polymorph phase diagram from first-principles quantum chemistry, mapping out the thermodynamic regions of stability for three polymorphs over the range 0–400 K and 0–6 GPa. The agreement between the predicted and experimental phase diagrams corresponds to predicting the relative polymorph free energies to within ∼0.5 kJ mol−1 accuracy, which is achieved by employing fragment-based second-order Møller–Plesset perturbation theory and coupled cluster theory plus a quasi-harmonic treatment of the phonons.
Whereas computational chemistry often strives for “chemical accuracy” of 1 kcal mol−1, organic polymorphs are often separated by 1–2 kJ mol−1 or less in the Gibbs free energies,15,16 and those free energies can exhibit considerable dependence on temperature and pressure. Successful organic polymorph phase diagram prediction requires sub-kJ mol−1 accuracy. Such accuracy has been achieved for the benzene lattice energy at 0 K,17 but accomplishing the same for the polymorph free energies across a range of temperatures and pressures is even harder.
While one can sometimes predict high-pressure phase transitions from enthalpy alone, predicting phase boundaries at low pressures requires modeling how the phonons and thermal expansion of the crystal impact the Gibbs free energy. First-principles phase diagram predictions are relatively common in inorganic materials,21,22 where rigid structures and dense crystal packing can simplify the free energy landscapes, but they are rarer for organic crystals. Most organic polymorph phase predictions have focused on transitions occurring at a specific temperature/pressure using a hamonic model,15,23 the quasi-harmonic approximation,16,24,25 or molecular dynamics.26–29 Complete phase diagram predictions have been most successful at high pressures30 or when using empirically fitted potentials.31,32
Methanol is one of the simplest organic molecules, but the occurrence of several phases over a relatively narrow range of temperatures and pressures, substantial variations in molar volume (∼22–36 cm3 mol−1) across these thermodynamic conditions, and strong cooperative hydrogen bonding effects make its phase diagram challenging to predict. Four crystalline phases of methanol have been reported (Fig. 1 and 2a). The orthorhombic α phase (P212121)33 occurs at low temperatures up to medium pressures. The disordered orthorhombic β phase (Cmc21)33 occurs at higher temperatures, while the triclinic γ phase (P)34 is stable at higher pressures and ambient temperatures. A high-pressure, low-temperature δ phase was recently reported,20 but its structure remains unknown.
Fig. 1 Experimental crystal structures for the α (left), disordered β (middle), and γ (right) polymorphs of methanol. |
Fig. 2 (a) Experimental phase diagram of methanol.18–20 See Section S1† for details. (b) Predicted phase diagram (solid lines) superimposed over the experimental one (dashed lines). |
The successful phase diagram prediction here is accomplished using fragment-based electronic structure methods coupled to a quasi-harmonic treatment35,36 of the Helmholtz vibrational free energy and thermal expansion. Specifically, the hybrid many-body interaction (HMBI) model partitions the crystal into molecules.37 The molecules and their short-range pairwise intermolecular interactions are computed quantum mechanically, while the long-range dimer and strong many-body interactions in methanol are approximated either with the AMOEBA polarizable force field38 or periodic Hartree–Fock (HF). Fragmenting allows the use of high-level electronic structure for the important monomer and dimer fragments, and lattice polarization effects are captured through the many-body treatment. The ability to converge crystal energetics with respect to the level of theory can be essential for resolving small polymorph energy differences.7,39 Quasi-harmonic fragment models can predict structural, thermochemical, and other molecular crystal properties in good agreement with experiment.25,40–43
The calculations reproduce various experimental structural and thermochemical data well, as shown in Table 1 and Sections S3 and S4.† For example, they predict an α-phase 0 K sublimation enthalpy of 45.0 kJ mol−1, versus 45.7 ± 0.3 kJ mol−1 experimentally.44 That value improves upon earlier predictions by up to several kJ mol−1.42,45,46 The predicted change in enthalpy associated with the α → β transition at 157 K and ambient pressure is 0.41 kJ mol−1, versus 0.64 kJ mol−1 experimentally.47 The experimentally observed decrease in the isobaric heat capacity across the α–β phase transition47 is reproduced qualitatively, though the heat capacities themselves are underestimated by 10–20% (Fig. 3). The underestimation of the heat capacity likely stems from a mixture of overestimation of the lattice mode phonon frequencies, the neglect of anharmonicity in the intramolecular modes (particularly methyl rotations), and underestimation of the thermal expansivity.42
Parameter | Experiment | Electronic CCSD(T)/CBS | Quasi-harmonic CCSD(T)/CBS |
---|---|---|---|
a METHOL04, ref. 33. b Ref. 44. c METHOL05, ref. 33. d METHOL03, ref. 34. | |||
α-Methanol (Z = 4, P2 1 2 1 2 1 , 122 K) | |||
V (Å3) | 207.04 | 191.80 | 208.26 |
a (Å) | 4.65 | 4.50 | 4.55 |
b (Å) | 4.93 | 4.80 | 5.00 |
c (Å) | 9.04 | 8.88 | 9.15 |
E lat (kJ mol−1) | — | 51.66 | 51.33 |
ΔHsub (0 K) (kJ mol−1) | 45.7 ± 0.3b | 44.70 | 44.95 |
β-Methanol (Z = 4, Cmc2 1 , 160 K) | |||
V (Å3) | 214.76 | 199.00 | 227.48 |
a (Å) | 6.40 | 6.06 | 6.68 |
b (Å) | 7.22 | 7.28 | 7.39 |
c (Å) | 4.65 | 4.51 | 4.61 |
E lat (kJ mol−1) | — | 50.82 | 50.54 |
ΔHsub (0 K) (kJ mol−1) | — | 44.43 | 44.71 |
γ-Methanol (Z = 6, P, 298 K, 4.0 GPa) | |||
V (Å3) | 236.80 | 235.18 | 230.90 |
a (Å) | 7.67 | 7.76 | 7.68 |
b (Å) | 4.41 | 4.34 | 4.33 |
c (Å) | 7.20 | 7.16 | 7.13 |
α (Å) | 88.10 | 85.82 | 85.78 |
β (Å) | 102.89 | 102.07 | 101.99 |
γ (Å) | 93.85 | 93.85 | 93.58 |
Fig. 3 Comparison between predicted and experimental47 isobaric heat capacities near the experimental α–β phase transition at 157 K. |
The phase diagram was constructed by computing G(T,p) for each phase over a range of temperatures and pressures and interpolating to find the phase equilibrium conditions. Fig. 4 plots sample free energy profiles at selected pressures. Fig. 2b overlays the predicted α, β, and γ coexistence curves onto the experimental phase diagram. Phase boundaries involving the unknown δ form were not computed. The predicted phase boundaries broadly reproduce the thermodynamic stability regions for the three polymorphs. The predictions underestimate the α–β transition temperature by ∼80 K near ambient pressure and ∼50 K or less at higher pressures. The α–γ transition pressure is predicted to within a fraction of a gigapascal. The predicted triple point for the three polymorphs occurs at 4.2 GPa and 250 K, versus 3.6 GPa and 290 K experimentally (errors of 0.6 GPa and 40 K).
For comparison, predictions of the methanol phase diagram with empirical force fields found no region of phase stability for the β phase whatsoever and predicted that the α–γ transition occurs above 10 GPa, instead of near 3.6 GPa.48,49 A carbon dioxide phase diagram predicted from density functional theory with the Perdew–Burke–Ernzerhof (PBE) functional underestimated the temperature-dependent phase transition in the 5–10 GPa regime by ∼200 K,30versus only ∼25 K errors for methanol in the same pressure regime here.
Further perspective on the quality of the predicted transition temperatures can be gleaned from sensitivity analysis (Fig. 5 and Section S5†). At ambient pressure, the α and β phases lie extremely close in free energy (ΔG = 0.24 kJ mol−1 at 0 K), and their free energies vary similarly with temperature (Fig. 4a). This translates to extreme sensitivity of the predicted phase transition temperature to small changes in the model. Increasing the predicted α–β phase transition temperature at ambient pressure from the current 81 K to the experimental 157 K would require stabilizing the α phase Gibbs free energy by a mere 0.4 kJ mol−1 relative to β. Similarly, a 0.5 kJ mol−1 change in the relative α–γ free energies would shift that coexistence curve to near-perfect agreement with experiment. In other words, the predicted phase boundaries here correspond to ∼0.5 kJ mol−1 accuracy or better in the relative free energies. At higher pressures, the increased importance of crystal packing density and diminished impact of thermal expansion make phase transitions easier to predict by reducing the sensitivity to small changes in the energies (Fig. 4c).
Fig. 5 Sensitivity of two phase coexistence curves to shifting the α-phase free energy by a small amount. |
Additional sensitivity analysis described in Section S5† reveals that the errors in the predicted transition temperatures stem largely from errors in the relative stabilities of the phases, while problems in the coexistence curve slopes arise from subtle variations in the curvatures of the electronic energy wells. The phase diagram is slightly less sensitive to errors in the phonons, though quasi-harmonic treatment of thermal expansion is key to reproducing the phase diagram even qualitatively.
The calculations here represent the disordered β phase via a single dynamically averaged structure. This is probably a reasonable approximation for the enthalpy, since experimental evidence33,50–56 largely points toward a thermally averaged structure and MP2 + AMOEBA calculations indicate that the disordered structures with out-of-plane distortion lie ∼1.5 kJ mol−1 higher than the averaged structure. On the other hand, the modeling approach may underestimate the entropic contribution from the disorder. For a static disorder with two possible orientations for each of the four molecules in the cell, one might simplistically estimate a configurational entropy stabilization for β methanol as Sconfig = Rln24 per cell, or 5.8 J mol−1 K−1 per molecule. That entropic stabilization would lower the α–β transition temperature by an additional ∼50 K at ambient pressure. Alternatively, the dynamically disordered scenario suggested by experiments and our calculations would manifest instead as a highly anharmonic motion that would not be captured by the quasi-harmonic phonon model here. See Section S4.2† for more details. A more careful treatment of the disorder28,57,58 would be a worthwhile future investigation.
Two other notable observations emerge from this study. First, it is common in crystal structure prediction to include Gibbs free energy estimates only in the final energy ranking steps. However, as shown in Fig. 6, the α phase exhibits a minimum on the free energy surface only below ∼240 K. Similarly, the γ phase only occurs at pressures above ∼1 GPa. These examples highlight how the potential and free energy landscapes can differ qualitatively, and that the thermodynamic conditions employed in crystal structure prediction can impact the structures that will be identified. Second, a candidate structure was recently proposed for the δ phase from crystal structure prediction.59 Calculations on this proposed structure predict its Gibbs free energy lies several kJ mol−1 above those of the α and γ phases throughout the relevant temperature and pressure ranges (Section S7†), suggesting that the proposed structure cannot account for the experimentally observed δ phase.
Discussing the future of chemistry in 1950, Robert Heinlein wrote “When chemistry becomes a discipline, mathematical chemists will design new materials, predict their properties, and tell engineers how to make them—without entering the laboratory. We've got a long way to go on that one.”61 Heinlein's third goal has proved particularly challenging for simulation, but the ability to predict a phase diagram ab initio represents a significant step toward the day when computation will be able to both predict what crystal structures might form and tell laboratory scientists how to make them.
As described in our recent study of α methanol,42 the crystal structures were relaxed using the hybrid many-body interaction (HMBI) fragment approach,37 employing density-fitted MP2 in the aug-cc-pVTZ basis set62 (Molpro 2012.1 (ref. 63)) for the one-body and short-range two-body terms and the AMOEBA force field38 (Tinker 6.2 (ref. 64)) for the long-range (>9–10 Å) pairwise and many-body contributions. Because the γ phase involves a larger unit cell (Z = 6) and low P symmetry, it was optimized in the smaller aug-cc-pVDZ basis set. Counterpoise corrections were applied for the dimer interactions throughout.
To evaluate the quasi-harmonic free energies, electronic energy vs. volume curves were mapped out for each polymorph, as shown in Fig. S2.† These curves were based on 15 structure optimizations each: the original optimization plus six more at positive pressures (compression branch) and eight more at negative pressures (expansion branch) using HMBI MP2/aug-cc-pVTZ + AMOEBA. The application of positive/negative pressure allows the crystals to expand or contract anisotropically.42 Single-point energies on these structures were then refined by replacing the AMOEBA many-body contributions with ones evaluated from periodic HF/pob-TZVP65 using CRYSTAL09 (ref. 66) (see Section S3†).
In addition, the impact of monomer and dimer correlation treatment at the MP2/aug-cc-pVTZ, MP2/CBS, and CCSD(T)/CBS levels of theory was considered (again using Molpro). Complete-basis set results were obtained from two-point aug-cc-pVTZ and aug-cc-pVQZ extrapolation.67,68
The resulting energy versus volume curves for each polymorph were fitted to a double Murnaghan equation of state—that is fitting the compression and expansion branches separately to the Murnaghan equation of state,
(1) |
Next, harmonic Γ-point phonons were computed at several different volumes at the HMBI MP2/aug-cc-pVTZ + AMOEBA level of theory (aug-cc-pVDZ for γ methanol). At each volume, the Helmholtz vibrational free energy is computed as a function of temperature according to
(2) |
The total Helmholtz free energy of the crystal is given as A(T,V) = E(V) + Avib(T,V). Using the thermodynamic relationship , the Gibbs free energy is evaluated as
(3) |
Examples of the resulting free energy curves are shown in Fig. 6. After mapping out the Gibbs free energies over a grid of temperatures and pressures, the coexistence curves (where ΔG = 0 between two phases) were identified via interpolation. See ref. 42 for further details.
Footnote |
† Electronic Supplementary Information (ESI) available: Discussion of the experimental uncertainties, additional modeling details, structure overlays, sensitivity analysis, and investigation of the δ phase structure. See DOI: 10.1039/c8sc01237g |
This journal is © The Royal Society of Chemistry 2018 |