Modelling temperature-dependent properties of polymorphic organic molecular crystals

First-Principles Models of Polymorphism of Pharmaceuticals: Maximizing the Accuracy-to-Cost Ratio

Accuracy and sophistication of in silico models of structure, internal dynamics, and cohesion of molecular materials at finite temperatures increase over time. Applicability limits of ab initio polymorph ranking that would be feasible at reasonable costs currently represent crystals of moderately sized molecules (less than 20 nonhydrogen atoms) and simple unit cells (containing rather only one symmetry-irreducible molecule). Extending the applicability range of the underlying first-principles methods to larger systems with a real-life significance, and enabling to perform such computations in a high-throughput regime represent additional challenges to be tackled in computational chemistry. This work presents a novel composite method that combines the computational efficiency of density-functional tight-binding (DFTB) methods with the accuracy of density-functional theory (DFT). Being rooted in the quasi-harmonic approximation, it uses a cheap method to perform all of the costly scans of how static and dynamic characteristics of the crystal vary with respect to its volume. Such data are subsequently corrected to agree with a higher-level model, which must be evaluated only at a single volume of the crystal. It thus enables predictions of structural, cohesive, and thermodynamic properties of complex molecular materials, such as pharmaceuticals or organic semiconductors, at a fraction of the original computational cost. As the composite model retains the solid physical background, it suffers from a minimum accuracy deterioration compared to the full treatment with the costly approach. The novel methodology is demonstrated to provide consistent results for the structural and thermodynamic properties of real-life molecular crystals and their polymorph ranking.

Thermodynamic Stability Is a Poor Indicator of Cocrystallization in Models of Organic Molecules

Cocrystallizing a given molecule with another can be useful for adjusting the physical properties of molecules in the solid state. However, most combinations of molecules do not readily cocrystallize but form either one-component crystals or amorphous solids. Computational methods of crystal structure prediction can, in principle, identify the thermodynamically stable cocrystal and thus predict if molecules will cocrystallize or not. However, the pronounced polymorphism and tendency of many organic molecules to form disordered solids suggest that kinetic factors can play an important role in cocrystallization. The question remains: if a binary system of molecules has a thermodynamically stable cocrystal, will it indeed cocrystallize? To address this question, we simulate the crystallization of more than 2600 distinct pairs of chiral model molecules of similar size in 2D and calculate accurate crystal energy landscapes for all of them. Our analysis shows that thermodynamic criteria alone are unreliable in the prediction of cocrystallization. While the vast majority of cocrystals that form in our simulations are thermodynamically favorable, most coformer systems that have a thermodynamically stable cocrystal do not cocrystallize. We furthermore show that cocrystallization rates increase 3-fold when coformers are used that do not form well-ordered single-component crystals. Our results suggest that kinetic factors of cocrystallization are decisive in many cases.

Frontiers of molecular crystal structure prediction for pharmaceuticals and functional organic materials

The reliability of organic molecular crystal structure prediction has improved tremendously in recent years. Crystal structure predictions for small, mostly rigid molecules are quickly becoming routine. Structure predictions for larger, highly flexible molecules are more challenging, but their crystal structures can also now be predicted with increasing rates of success. These advances are ushering in a new era where crystal structure prediction drives the experimental discovery of new solid forms. After briefly discussing the computational methods that enable successful crystal structure prediction, this perspective presents case studies from the literature that demonstrate how state-of-the-art crystal structure prediction can transform how scientists approach problems involving the organic solid state. Applications to pharmaceuticals, porous organic materials, photomechanical crystals, organic semi-conductors, and nuclear magnetic resonance crystallography are included. Finally, efforts to improve our understanding of which predicted crystal structures can actually be produced experimentally and other outstanding challenges are discussed.

Polymorph screening at surfaces of a benzothienobenzothiophene derivative: discovering new solvate forms

The discovery of new polymorphs opens up unique applications for molecular materials since their physical properties are predominantly influenced by the crystal structure type. The deposition of molecules at surfaces offers great potential in the variation of the crystallization conditions, thereby allowing access to unknown polymorphs. With our surface crystallization approach, four new phases are found for an oligoethylene glycol-benzothienobenzothiophene molecule, and none of these phases could be identified via classical polymorph screening. The corresponding crystal lattices of three of the new phases were obtained via X-ray diffraction (XRD). Based on the volumetric considerations together with X-ray fluorescence and Raman spectroscopy data, the phases are identified as solvates containing one, two or three solvent molecules per molecule. The strong interaction of dichloromethane with the oligoethylene glycol side chains of the molecules may be responsible for the formation of the solvates. Temperature-dependent XRD reveals the low thermal stability of the new phases, contrary to the thermodynamically stable bulk form. Nevertheless, the four solvates are stable under ambient conditions for at least two years. This work illustrates that defined crystallization at surfaces enables access to multiple solvates of a given material through precise and controlled variations in the crystallization kinetics.

A general force field by machine learning on experimental crystal structures. Calculations of intermolecular Gibbs energy with FlexCryst

Machine learning was employed on the experimental crystal structures of the Cambridge Structural Database (CSD) to derive an intermolecular force field for all available types of atoms (general force field). The obtained pairwise interatomic potentials of the general force field allow for the fast and accurate calculation of intermolecular Gibbs energy. The approach is based on three postulates regarding Gibbs energy: the lattice energy must be below zero, the crystal structure must be a local minimum, and, if available, the experimental and the calculated lattice energy must coincide. The parametrized general force field was then validated regarding these three conditions. First, the experimental lattice energy was compared with the calculated energies. The observed errors were found to be in the order of experimental errors. Second, Gibbs lattice energy was calculated for all structures available in the CSD. Their energy values were found to be below zero in 99.86% of the cases. Finally, 500 random structures were minimized, and the change in density and energy was examined. The mean error in the case of density was below 4.06%, and for energy it was below 5.7%. The obtained general force field calculated Gibbs lattice energies of 259 041 known crystal structures within a few hours. Since Gibbs energy defines the reaction energy, the calculated energy can be used to predict chemical-physical properties of crystals, for instance, the formation of co-crystals, polymorph stability and solubility.

Discovery of a new polymorph of clotrimazole through melt crystallization: Understanding nucleation and growth kinetics

Clotrimazole (CMZ) is a classical antifungal drug for studying crystallization. In this study, a new CMZ polymorph (Form 2) was discovered during the process of nucleation and growth rate determination in the melt. High-quality single crystals were grown from melt microdroplets to determine the crystal structure by x-ray diffraction. Form 2 is metastable and exhibits a disordered structure. The crystal nucleation and growth kinetics of the two CMZ polymorphs were systematically measured. Form 2 nucleates and grows faster than the existing form (Form 1). The maximum nucleation rate of Forms 1 and 2 was observed at 50 °C (1.07 T_g). The summary of the maximum nucleation rate temperature of CMZ and the other six organic compounds indicates that nucleation near T_g in the supercooled liquid is a useful approach to discovering new polymorphs. This study is relevant for the discovering new drug polymorphs through an understanding of nucleation and growth kinetics during melt crystallization.

Application of gas chromatographic data and 2D molecular descriptors for accurate global mobility potential prediction

Mobility is a key feature affecting the environmental fate, which is of particular importance in the case of persistent organic pollutants (POPs) and emerging pollutants (EPs). In this study, the global mobility classification artificial neural networks-based models employing GC retention times (RT) and 2D molecular descriptors were constructed and validated. The high usability of RT was confirmed based on the feature selection step performed using the multivariate adaptive regression splines (MARS) tool. Although RT was found to be the most important, according to Kruskal-Wallis ANOVA analysis, it is insufficient to build a robust model, which justifies the need to expand the input layer with 2D descriptors. Therefore the following molecular descriptors: MPC10, WTPT-2, AATS8s, minaaCH, GATS7c, RotBtFrac, ATSC7v and ATSC1p, which were characterized by a high predicting potential were used to improve the classification performance. As a result of machine learning procedure ten of the most accurate neural networks were selected. The external validation showed that the final models are characterized by a high general accuracy score (85.71-96.43%). The high predicting abilities were also confirmed by the micro-averaged Matthews correlation coefficient (MAMCC) (0.73-0.88). To evaluate the applicability of the models, new retention times of selected POPs and EPs including pesticides, polycyclic aromatic hydrocarbons, pharmaceuticals, fragrances and personal care products were measured and used for mobility prediction. Further, the classifiers were used for photodegradation and chlorination products of two popular sunscreen agents, 2-ethyl-hexyl-4-methoxycinnamate and 2-ethylhexyl 4-(dimethylamino)benzoate.

The antisolvent coprecipitation method for enhanced bioavailability of poorly water-soluble drugs

In recent years, poorly water-soluble drug candidates in the drug development pipeline have been a challenging issue for the pharmaceutical industry. Many delivery systems such as nanocrystals, cocrystals, nanoparticles, and amorphous solid dispersions (ASDs) have been developed to overcome these problems. A large number of methods are utilized to realize the above delivery systems. Among all the preparation methods, the antisolvent coprecipitation method is a relatively simple, cost-effective method, offering many advantages over conventional methods. An overview of recent developments for each solubility enhancement approach using the antisolvent coprecipitation method is presented. This current review details a comprehensive overview of the antisolvent coprecipitation process and its properties, as well as the fundamentals for enhancing the solubility and bioavailability of poorly water-soluble drugs by nanotization, polymorph control with polymers and/or surfactants. Furthermore, this review also presents insights into the factors affecting the antisolvent coprecipitation process.

Performance of Dispersion-Inclusive Density Functional Theory Methods for Energetic Materials

Molecular crystals of energetic materials (EMs) are denser than typical molecular crystals and are characterized by distinct intermolecular interactions between nitrogen-containing moieties. To assess the performance of dispersion-inclusive density functional theory (DFT) methods, we have compiled a data set of experimental sublimation enthalpies of 31 energetic materials. We evaluate the performance of three methods: the semilocal Perdew-Burke-Ernzerhof (PBE) functional coupled with the pairwise Tkatchenko-Scheffler (TS) dispersion correction, PBE with the many-body dispersion (MBD) method, and the PBE-based hybrid functional (PBE0) with MBD. Zero-point energy contributions and thermal effects are described using the quasi-harmonic approximation (QHA), including explicit treatment of thermal expansion, which we find to be non-negligible for EMs. The lattice energies obtained with PBE0+MBD are the closest to experimental sublimation enthalpies with a mean absolute error of 9.89 kJ/mol. However, the state-of-the-art treatment of vibrational and thermal contributions makes the agreement with experiment worse. Pressure-volume curves are also examined for six representative materials. For pressure-volume curves, all three methods provide reasonable agreement with experimental data with mean absolute relative errors of 3% or less. Most of the intermolecular interactions typical of EMs, namely nitro-amine, nitro-nitro, and nitro-hydrogen interactions, are more sensitive to the choice of the dispersion method than to the choice of the exchange-correlation functional. The exception is π-π stacking interactions, which are also very sensitive to the choice of the functional. Overall, we find that PBE+TS, PBE+MBD, and PBE0+MBD do not perform as well for energetic materials as previously reported for other classes of molecular crystals. This highlights the importance of testing dispersion-inclusive DFT methods for diverse classes of materials and the need for further method development.

Polymorphism in Ionic Cocrystals Comprising Lithium Salts and l-Proline

The occurrence of polymorphism in ionic cocrystals formed by two lithium salts, lithium salicylate (LIS) and lithium 4-methoxybenzoate (L4M), and l-proline (PRO) has been investigated. The previously reported monoclinic form of the 1:1 cocrystal of LIS and PRO, LISPRO(α), and a new thermodynamically stable orthorhombic polymorph, LISPRO(β), were prepared and characterized. The two polymorphs form square grid, sql, topology coordination networks and differ mainly in the conformation of the salicylate ions and positioning of the sql nets. LISPRO(α) was observed to transform to LISPRO(β) under slurry conditions. The 1:1 ionic cocrystal of L4M and PRO (L4MPRO) was found to form three polymorphs. Apart from the previously reported orthorhombic crystal form, L4MPRO(α), two new monoclinic crystal forms, L4MPRO(β) and L4MPRO(γ), were obtained by modifying crystallization conditions. The new polymorphs were found to be metastable, undergoing transformations to L4MPRO(α) upon exposure to humidity. Experimental conditions that induce transformations between the polymorphs of LISPRO and L4MPRO are detailed, and the structural differences between the polymorphs are discussed in the broader context of polymorphism.

Thermal Properties and Ionic Conductivity of Tetra-n-Butylammonium Perchlorate

The thermal parameters of the phase transitions and transport properties of tetra-n-butylammonium (TBA) perchlorate (n-C4H9)4NClO4 (TBAClO4) were investigated. TBAClO4 has a polymorphous transition at 330 K and melts at 487 K. The structure of the high-temperature (HT) phase belongs to cubic symmetry and is similar to the HT phases of TBABF4 and TBAI salts. The conductivity parameters of the low-temperature and HT phases of TBAClO4 were determined from the Arrhenius plots. The thermodynamic parameters and transport properties of TBAClO4 were compared with those of other TBA salts having isostructural HT phases. The polymorphous phase transition entropy was found to correlate with the conductivity of HT phases of TBA salts; TBAClO4 has the lowest conductivity compared to TBABF4 and TBAI salts.

A density-functional benchmark of vibrational free-energy corrections for molecular crystal polymorphism

Many crystal structure prediction protocols only concern themselves with the electronic energy of molecular crystals. However, vibrational contributions to the free energy ( Fvib) can be significant in determining accurate stability rankings for crystal candidates. While force-field studies have been conducted to gauge the magnitude of these free-energy corrections, highly accurate results from quantum mechanical methods, such as density-functional theory (DFT), are desirable. Here, we introduce the PV17 set of 17 polymorphic pairs of organic molecular crystals, for which plane wave DFT is used to calculate the vibrational free energies and free-energy differences (Δ Fvib) between each pair. Our DFT results confirm that the vibrational free-energy corrections are small, having a mean value of 1.0 kJ/mol and a maximum value of 2.3 kJ/mol for the PV17 set. Furthermore, we assess the accuracy of a series of lower-cost DFT, semi-empirical, and force-field models for computing Δ Fvib that have been proposed in the literature. It is found that calculating Fvib using the Γ-point frequencies does not provide Δ Fvib values of sufficiently high quality. In addition, Δ Fvib values calculated using various approximate methods have mean absolute errors relative to our converged DFT results of equivalent or larger magnitude than the vibrational free-energy corrections themselves. Thus, we conclude that, in a crystal structure prediction protocol, it is preferable to forego the inclusion of vibrational free-energy corrections than to estimate them with any of the approximate methods considered here.

Emerging Landscape of Computational Modeling in Pharmaceutical Development

Computational chemistry applications have become an integral part of the drug discovery workflow over the past 35 years. However, computational modeling in support of drug development has remained a relatively uncharted territory for a significant part of both academic and industrial communities. This review considers the computational modeling workflows for three key components of drug preclinical and clinical development, namely, process chemistry, analytical research and development, as well as drug product and formulation development. An overview of the computational support for each step of the respective workflows is presented. Additionally, in context of solid form design, special consideration is given to modern physics-based virtual screening methods. This covers rational approaches to polymorph, coformer, counterion, and solvent virtual screening in support of solid form selection and design.

Crystal Structure Prediction Methods for Organic Molecules: State of the Art

The prediction of the crystal structures that a given organic molecule is likely to form is an important theoretical problem of significant interest for the pharmaceutical and agrochemical industries, among others. As evidenced by a series of six blind tests organized over the past 2 decades, methodologies for crystal structure prediction (CSP) have witnessed substantial progress and have now reached a stage of development where they can begin to be applied to systems of practical significance. This article reviews the state of the art in general-purpose methodologies for CSP, placing them within a common framework that highlights both their similarities and their differences. The review discusses specific areas that constitute the main focus of current research efforts toward improving the reliability and widening applicability of these methodologies, and offers some perspectives for the evolution of this technology over the next decade.

Exploration and Optimization in Crystal Structure Prediction: Combining Basin Hopping with Quasi-Random Sampling

We describe the implementation of a Monte Carlo basin hopping (BH) global optimization procedure for the prediction of molecular crystal structures. The BH method is combined with quasi-random (QR) structure generation in a hybrid method for crystal structure prediction, QR-BH, which combines the low-discrepancy sampling provided by QR sequences with BH efficiency at locating low energy structures. Through tests on a set of single-component molecular crystals and co-crystals, we demonstrate that QR-BH provides faster location of low energy structures than pure QR sampling, while maintaining the efficient location of higher energy structures that are important for identifying important polymorphs.

The unexpected discovery of the ninth polymorph of tolfenamic acid

A new polymorph of tolfenamic acid, form IX, has been crystallised from a simple cooling crystallisation experiment raising the question as to why this polymorph had never been reported before.

Improved quantitative crystal-structure comparison using powder diffractograms via anisotropic volume correction

A new anisotropic volume correction improves quantitative crystal structure comparison. Benchmarking against the 6th crystal structure prediction blind test data results in identification of two previously uncredited matching structures.

Suppression of isotopic polymorphism

Crystallisation at pressure overcomes the effect of isotopic polymorphism in the methylpyridine pentachlorophenol co-crystal. Though the hydrogenated Cc polymorph can only be obtained at pressure, it is stable on recovery to ambient conditions.

Diabat method for polymorph free energies: Extension to molecular crystals

Lattice-switch Monte Carlo and the related diabat methods have emerged as efficient and accurate ways to compute free energy differences between polymorphs. In this work, we introduce a one-to-one mapping from the reference positions and displacements in one molecular crystal to the positions and displacements in another. Two features of the mapping facilitate lattice-switch Monte Carlo and related diabat methods for computing polymorph free energy differences. First, the mapping is unitary so that its Jacobian does not complicate the free energy calculations. Second, the mapping is easily implemented for molecular crystals of arbitrary complexity. We demonstrate the mapping by computing free energy differences between polymorphs of benzene and carbamazepine. Free energy calculations for thermodynamic cycles, each involving three independently computed polymorph free energy differences, all return to the starting free energy with a high degree of precision. The calculations thus provide a force field independent validation of the method and allow us to estimate the precision of the individual free energy differences.

Reduced-cost supercell approach for computing accurate phonon density of states in organic crystals

Phonon contributions to organic crystal structures and thermochemical properties can be significant, but computing a well-converged phonon density of states with lattice dynamics and periodic density functional theory (DFT) is often computationally expensive due to the need for large supercells. Using semi-empirical methods like density functional tight binding (DFTB) instead of DFT can reduce the computational costs dramatically, albeit with noticeable reductions in accuracy. This work proposes approximating the phonon density of states via a relatively inexpensive DFTB supercell treatment of the phonon dispersion that is then corrected by shifting the individual phonon modes according to the difference between the DFT and DFTB phonon frequencies at the Γ-point. The acoustic modes are then computed at the DFT level from the elastic constants. In several small-molecule crystal test cases, this combined approach reproduces DFT thermochemistry with kJ/mol accuracy and 1-2 orders of magnitude less computational effort. Finally, this approach is applied to computing the free energy differences between the five crystal polymorphs of oxalyl dihydrazide.

Structure prediction of crystals, surfaces and nanoparticles

We review the current techniques used in the prediction of crystal structures and their surfaces and of the structures of nanoparticles. The main classes of search algorithm and energy function are summarized, and we discuss the growing role of methods based on machine learning. We illustrate the current status of the field with examples taken from metallic, inorganic and organic systems. This article is part of a discussion meeting issue 'Dynamic in situ microscopy relating structure and function'.

Statistical Mechanical Approximations to More Efficiently Determine Polymorph Free Energy Differences for Small Organic Molecules

Methods to efficiently determine the relative stability of polymorphs of organic crystals are highly desired in crystal structure predictions (CSPs). Current methodologies include calculating the free energy of static lattice phonons, quasi-harmonic approximations (QHA), and computing the full thermodynamic cycle using replica exchange molecular dynamics (REMD). We found that 13 out of the 29 systems minimized from experimental crystal structures restructured to a lower energy minimum when heated and annealed using REMD, a phenomenon that QHA alone cannot capture. Here, we present a series of methods that are intermediate in accuracy and expense between QHA and computing the full thermodynamic cycle, which can save 42-80% of the computational cost and introduces, on this benchmark, a relatively small (0.16 ± 0.04 kcal/mol) error relative to the full thermodynamic cycle. In particular, a method that Boltzmann weights harmonic free energies from along the trajectory of REMD replica exchange appears to be an appropriate intermediate between QHA and the full thermodynamic cycle using MD when screening crystal polymorph stability.

Pharmaceutical Hydrates Analysis-Overview of Methods and Recent Advances

This review discusses a set of instrumental and computational methods that are used to characterize hydrated forms of APIs (active pharmaceutical ingredients). The focus has been put on highlighting advantages as well as on presenting some limitations of the selected analytical approaches. This has been performed in order to facilitate the choice of an appropriate method depending on the type of the structural feature that is to be analyzed, that is, degree of hydration, crystal structure and dynamics, and (de)hydration kinetics. The presented techniques include X-ray diffraction (single crystal X-ray diffraction (SCXRD), powder X-ray diffraction (PXRD)), spectroscopic (solid state nuclear magnetic resonance spectroscopy (ssNMR), Fourier-transformed infrared spectroscopy (FT-IR), Raman spectroscopy), thermal (differential scanning calorimetry (DSC), thermogravimetric analysis (TGA)), gravimetric (dynamic vapour sorption (DVS)), and computational (molecular mechanics (MM), Quantum Mechanics (QM), molecular dynamics (MD)) methods. Further, the successful applications of the presented methods in the studies of hydrated APIs as well as studies on the excipients' influence on these processes have been described in many examples.

The effective volumes of waters of crystallization: general organic solids

Using a list of compatible hydrate/anhydrate pairs prepared by van de Streek and Motherwell [CrystEngComm (2007), 9, 55-64], we have examined the effective volume per water of crystallization for 179 pairs of organic solids using current data from the Cambridge Crystallographic Structural Database (CSD). The effective volume is the difference per water molecule between the asymmetric unit volumes of the hydrate and parent anhydrate, and has the mean value 24 ų. The conformational changes in the reference molecule between the hydrate and its anhydrate are shown in two figures: one for a relatively rigid standard organic molecule and (in the supplementary file) one for a more flexible linear molecule. Using data from Nyman and Day [Phys. Chem. Chem. Phys. (2016), 18, 31132-31143], we have also established a generic volumetric coefficient of thermal expansion of organic solids with a value of 147 ± 56 × 10^-6 K^-1. There is a significant number of outliers to the data, negative, near zero, and large and positive. Some explanation for the existence of these outliers is attempted.

Sublimation Properties of α,ω-Diamines Revisited from First-Principles Calculations

Sublimation enthalpies of alkane-α,ω-diamines exhibit an odd-even pattern within their homologous series. First-principles calculations coupled with the quasi-harmonic approximation for crystals and with the conformation mixing model for the ideal gas are used to explain this phenomenon from the theoretical point of view. Crystals of the odd and even alkane-α,ω-diamines distinctly differ in their packing motifs. However, first-principles calculations indicate that it is a delicate interplay of the cohesive forces, phonons, molecular vibrations and conformational equilibrium which governs the odd-even pattern of the sublimation enthalpies within the homologous series. High molecular flexibility of the alkane-α,ω-diamines predetermines higher sensitivity of the computational model to the quality of the optimized geometries and relative conformational energies. Performance of high-throughput computational methods, such as the density functional tight binding (DFTB, GFN2-xTB) and the explicitly correlated dispersion-corrected Møller-Plesset perturbative method (MP2C-F12), are benchmarked against the consistent state-of-the-art calculations of conformational energies and interaction energies, respectively.

Evolutionary chemical space exploration for functional materials: computational organic semiconductor discovery

Computational methods, including crystal structure and property prediction, have the potential to accelerate the materials discovery process by enabling structure prediction and screening of possible molecular building blocks prior to their synthesis. However, the discovery of new functional molecular materials is still limited by the need to identify promising molecules from a vast chemical space. We describe an evolutionary method which explores a user specified region of chemical space to identify promising molecules, which are subsequently evaluated using crystal structure prediction. We demonstrate the methods for the exploration of aza-substituted pentacenes with the aim of finding small molecule organic semiconductors with high charge carrier mobilities, where the space of possible substitution patterns is too large to exhaustively search using a high throughput approach. The method efficiently explores this large space, typically requiring calculations on only ∼1% of molecules during a search. The results reveal two promising structural motifs: aza-substituted naphtho[1,2-a]anthracenes with reorganisation energies as low as pentacene and a series of pyridazine-based molecules having both low reorganisation energies and high electron affinities.

Overcoming the difficulties of predicting conformational polymorph energetics in molecular crystals via correlated wavefunction methods

Molecular crystal structure prediction is increasingly being applied to study the solid form landscapes of larger, more flexible pharmaceutical molecules. Despite many successes in crystal structure prediction, van der Waals-inclusive density functional theory (DFT) methods exhibit serious failures predicting the polymorph stabilities for a number of systems exhibiting conformational polymorphism, where changes in intramolecular conformation lead to different intermolecular crystal packings. Here, the stabilities of the conformational polymorphs of o-acetamidobenzamide, ROY, and oxalyl dihydrazide are examined in detail. DFT functionals that have previously been very successful in crystal structure prediction perform poorly in all three systems, due primarily to the poor intramolecular conformational energies, but also due to the intermolecular description in oxalyl dihydrazide. In all three cases, a fragment-based dispersion-corrected second-order Møller-Plesset perturbation theory (MP2D) treatment of the crystals overcomes these difficulties and predicts conformational polymorph stabilities in good agreement with experiment. These results highlight the need for methods which go beyond current-generation DFT functionals to make crystal polymorph stability predictions truly reliable.

Intermolecular Interactions in Molecular Organic Crystals upon Relaxation of Lattice Parameters

Crystal structure prediction is based on the assumption that the most thermodynamically stable structure will crystallize first. The existence of other structures such as polymorphs or from counterenantiomers requires an accurate calculation of the electronic energy. Using atom-centered Gaussian basis functions in periodic Density Functional Theory (DFT) calculations in Turbomole, the performance of two dispersion-corrected functionals, PBE-D3 and B97-D, is assessed for molecular organic crystals of the X23 benchmark set. B97-D shows a MAE (mean absolute error) of 4 kJ/mol, compared to 9 kJ/mol for PBE-D3. A strategy for the convergence of lattice energies towards the basis set limit is outlined. A simultaneous minimization of molecular structures and lattice parameters shows that both methods are able to reproduce experimental unit cell parameters to within 4–5%. Calculated lattice energies, however, deviate slightly more from the experiment, i.e., by 0.4 kJ/mol after unit cell optimization for PBE-D3 and 0.5 kJ/mol for B97-D. The accuracy of the calculated lattice energies compared to the experimental values demonstrates the ability of current DFT methods to assist in the quest for possible polymorphs and enantioselective crystallization processes.

Improving Predicted Nuclear Magnetic Resonance Chemical Shifts Using the Quasi-Harmonic Approximation

Ab initio nuclear magnetic resonance chemical shift prediction plays an important role in the determination or validation of crystal structures. The ability to predict chemical shifts more accurately can translate to increased confidence in the resulting chemical shift or structural assignments. Standard electronic structure predictions for molecular crystal structures neglect thermal expansion, which can lead to an appreciable underestimation of the molar volumes. This study examines this volume error and its impact on 68 ¹³C- and 28 ¹⁵N-predicted chemical shifts taken from 20 molecular crystals. It assesses the ability to recover more realistic room-temperature crystal structures using the quasi-harmonic approximation and how refining the structures impacts the chemical shifts. Several pharmaceutical molecular crystals are also examined in more detail. On the whole, accounting for quasi-harmonic expansion changes the ¹³C and ¹⁵N chemical shifts by 0.5 and 1.0 ppm on average. This, in turn, reduces the root-mean-square errors relative to experiment by 0.3 ppm for ¹³C and 0.7 ppm for ¹⁵N. Although the statistical impacts are modest, changes in individual chemical shifts can reach multiple ppm. Accounting for thermal expansion in molecular crystal chemical shift prediction may not be needed routinely, but the systematic trend toward improved accuracy with the experiment could be useful in cases where discrimination between structural candidates is challenging, as in the pharmaceutical theophylline.

X-ray diffraction data as a source of the vibrational free-energy contribution in polymorphic systems

In this contribution we attempt to answer a general question: can X-ray diffraction data combined with theoretical computations be a source of information about the thermodynamic properties of a given system? Newly collected sets of high-quality multi-temperature single-crystal X-ray diffraction data and complementary periodic DFT calculations of vibrational frequencies and normal mode vectors at the Γ point on the yellow and white polymorphs of di-methyl 3,6-di-chloro-2,5-di-hydroxy-terephthalate are combined using two different approaches, aiming to obtain thermodynamic properties for the two compounds. The first approach uses low-frequency normal modes extracted from multi-temperature X-ray diffraction data (normal coordinate analysis), while the other uses DFT-calculated low-frequency normal mode in the refinement of the same data (normal mode refinement). Thermodynamic data from the literature [Yang et al. (1989), Acta Cryst. B45, 312-323] and new periodic ab initio DFT supercell calculations are used as a reference point. Both approaches tested in this work capture the most essential features of the systems: the polymorphs are enantiotropically related, with the yellow form being the thermodynamically stable system at low temperature, and the white form at higher temperatures. However, the inferred phase transition temperature varies between different approaches. Thanks to the application of unconventional methods of X-ray data refinement and analysis, it was additionally found that, in the case of the yellow polymorph, anharmonicity is an important issue. By discussing contributions from low- and high-frequency modes to the vibrational entropy and enthalpy, the importance of high-frequency modes is highlighted. The analysis shows that larger anisotropic displacement parameters are not always related to the polymorph with the higher vibrational entropy contribution.

Revised values for the X23 benchmark set of molecular crystals

A revised reference value set for molecular crystals: X23b; new cell volumes and lattice energies including volumetric expansion due to zero-point energy and thermal effects.

Accuracy and reproducibility in crystal structure prediction: the curious case of ROY

Because of excessive electron delocalization, the polymorphs of ROY constitute a surprisingly challenging system for crystal structure prediction.

Thermal Gradient Approach for the Quasi-harmonic Approximation and Its Application to Improved Treatment of Anisotropic Expansion

We present a novel approach to efficiently implement thermal expansion in the quasi-harmonic approximation (QHA) for both isotropic and more importantly, anisotropic expansion. In this approach, we rapidly determine a crystal's equilibrium volume and shape at a given temperature by integrating along the gradient of expansion from 0 Kelvin up to the desired temperature. We compare our approach to previous isotropic methods that rely on a brute-force grid search to determine the free energy minimum, which is infeasible to carry out for anisotropic expansion, as well as quasi-anisotropic approaches that take into account the contributions to anisotropic expansion from the lattice energy. We compare these methods for experimentally known polymorphs of piracetam and resorcinol and show that both isotropic methods agree to within error up to 300 K. Using the Grüneisen parameter causes up to 0.04 kcal/mol deviation in the Gibbs free energy, but for polymorph free energy differences there is a cancellation in error with all isotropic methods within 0.025 kcal/mol at 300 K. Anisotropic expansion allows the crystals to relax into lattice geometries 0.01-0.23 kcal/mol lower in energy at 300 K relative to isotropic expansion. For polymorph free energy differences all QHA methods produced results within 0.02 kcal/mol of each other for resorcinol and 0.12 kcal/mol for piracetam, the two molecules tested here, demonstrating a cancellation of error for isotropic methods. We also find that with expansion in more than a single volume variable, there is a non-negligible rate of failure of the basic approximations of QHA. Specifically, while expanding into new harmonic modes as the box vectors are increased, the system often falls into alternate, structurally distinct harmonic modes unrelated by continuous deformation from the original harmonic mode.

A thermodynamically consistent phase diagram of a trimorphic pharmaceutical, l-tyrosine ethyl ester, based on limited experimental data

Crystalline polymorphs possess different physical properties, and phase changes between those polymorphs may affect the properties of engineered materials such as drugs. This is very well illustrated by the large effort that is put into the capability to predict phase behaviour of pharmaceuticals to avoid the unexpected appearance of different crystal forms. Much progress has been made, but one of the remaining challenges is (the accuracy in) the prediction of phase behaviour as a function of temperature. Obviously, predictions should at a certain point be verified against experimental data; however, it may not always be easy to elucidate the phase behaviour of a given compound experimentally, because thermodynamically and kinetically controlled phenomena occur in a convoluted fashion in experimental data. The present paper discusses the trimorphism of l-tyrosine ethyl ester as an example case of how experimental data in combination with the thermodynamic tenets lead to a consistent phase diagram, which can be used as the basis for pharmaceutical formulations and for comparison with polymorph predictions by computer. The positions of the two-phase equilibria I-II, I-III, and I-L have been obtained experimentally. Using the Clapeyron equation and the alternation rule, it has been shown how the positions of the other equilibria II-L, III-L, and II-III can be deduced in combination with the stability rankings of the phases and the phase equilibria. The experimental data have been obtained by synchrotron X-ray diffraction, Raman spectroscopy, and thermal analysis as a function of pressure and temperature. Furthermore, laboratory X-ray diffraction as a function of temperature and differential scanning calorimetry have been used. At room temperature, form II is the most stable phase, which remains stable with increasing pressure, as it possesses the smallest specific volume. Form I becomes stable above 33 °C (306 K), but with increasing pressure it turns into form III. On thermodynamic grounds, form III is expected to have a stable domain at very low temperatures.

Control and prediction of the organic solid state: a challenge to theory and experiment†

The ability of theoretical chemists to quantitatively model the weak forces between organic molecules is being exploited to predict their crystal structures and estimate their physical properties. Evolving crystal structure prediction methods are increasingly being used to aid the design of organic functional materials and provide information about thermodynamically plausible polymorphs of speciality organic materials to aid, for example, pharmaceutical development. However, the increasingly sophisticated experimental studies for detecting the range of organic solid-state behaviours provide many challenges for improving quantitative theories that form the basis for the computer modelling. It is challenging to calculate the relative thermodynamic stability of different organic crystal structures, let alone understand the kinetic effects that determine which polymorphs can be observed and are practically important. However, collaborations between experiment and theory are reaching the stage of devising experiments to target the first crystallization of new polymorphs or create novel organic molecular materials.

Ab initio prediction of the polymorph phase diagram for crystalline methanol

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.

Evaluating the Energetic Driving Force for Cocrystal Formation

We present a periodic density functional theory study of the stability of 350 organic cocrystals relative to their pure single-component structures, the largest study of cocrystals yet performed with high-level computational methods. Our calculations demonstrate that cocrystals are on average 8 kJ mol^-1 more stable than their constituent single-component structures and are very rarely (<5% of cases) less stable; cocrystallization is almost always a thermodynamically favorable process. We consider the variation in stability between different categories of systems-hydrogen-bonded, halogen-bonded, and weakly bound cocrystals-finding that, contrary to chemical intuition, the presence of hydrogen or halogen bond interactions is not necessarily a good predictor of stability. Finally, we investigate the correlation of the relative stability with simple chemical descriptors: changes in packing efficiency and hydrogen bond strength. We find some broad qualitative agreement with chemical intuition-more densely packed cocrystals with stronger hydrogen bonding tend to be more stable-but the relationship is weak, suggesting that such simple descriptors do not capture the complex balance of interactions driving cocrystallization. Our conclusions suggest that while cocrystallization is often a thermodynamically favorable process, it remains difficult to formulate general rules to guide synthesis, highlighting the continued importance of high-level computation in predicting and rationalizing such systems.

The effect of amino acid backbone length on molecular packing: crystalline tartrates of glycine, β-alanine, γ-aminobutyric acid (GABA) and DL-α-aminobutyric acid (AABA)

We report a novel 1:1 cocrystal of β-alanine with DL-tartaric acid, C₃H₇NO₂·C₄H₆O₆, (II), and three new molecular salts of DL-tartaric acid with β-alanine {3-azaniumylpropanoic acid-3-azaniumylpropanoate DL-tartaric acid-DL-tartrate, [H(C₃H₇NO₂)₂]⁺·[H(C₄H₅O₆)₂]^-, (III)}, γ-aminobutyric acid [3-carboxypropanaminium DL-tartrate, C₄H₁₀NO₂⁺·C₄H₅O₆^-, (IV)] and DL-α-aminobutyric acid {DL-2-azaniumylbutanoic acid-DL-2-azaniumylbutanoate DL-tartaric acid-DL-tartrate, [H(C₄H₉NO₂)₂]⁺·[H(C₄H₅O₆)₂]^-, (V)}. The crystal structures of binary crystals of DL-tartaric acid with glycine, (I), β-alanine, (II) and (III), GABA, (IV), and DL-AABA, (V), have similar molecular packing and crystallographic motifs. The shortest amino acid (i.e. glycine) forms a cocrystal, (I), with DL-tartaric acid, whereas the larger amino acids form molecular salts, viz. (IV) and (V). β-Alanine is the only amino acid capable of forming both a cocrystal [i.e. (II)] and a molecular salt [i.e. (III)] with DL-tartaric acid. The cocrystals of glycine and β-alanine with DL-tartaric acid, i.e. (I) and (II), respectively, contain chains of amino acid zwitterions, similar to the structure of pure glycine. In the structures of the molecular salts of amino acids, the amino acid cations form isolated dimers [of β-alanine in (III), GABA in (IV) and DL-AABA in (V)], which are linked by strong O-H...O hydrogen bonds. Moreover, the three crystal structures comprise different types of dimeric cations, i.e. (A...A)⁺ in (III) and (V), and A⁺...A⁺ in (IV). Molecular salts (IV) and (V) are the first examples of molecular salts of GABA and DL-AABA that contain dimers of amino acid cations. The geometry of each investigated amino acid (except DL-AABA) correlates with the melting point of its mixed crystal.

Crystal structure prediction is changing from basic science to applied technology

Prediction of true polymorphs as dynamic ensembles in contrast to hypothetical static crystal structures.

Crystal structure prediction of flexible pharmaceutical-like molecules: density functional tight-binding as an intermediate optimisation method and for free energy estimation

Successful methodologies for theoretical crystal structure prediction (CSP) on flexible pharmaceutical-like organic molecules explore the lattice energy surface to find a set of plausible crystal structures. The initial search stages of CSP studies use relatively simple lattice energy approximations as hundreds of thousands of minima have to be considered. These generated crystal structures often have poor molecular geometries, as well as inaccurate lattice energy rankings, and performing reasonably accurate but computationally affordable optimisations of the crystal structures generated in a search would be highly desirable. Here, we seek to explore whether semi-empirical quantum-mechanical methods can perform this task. We employed the dispersion-corrected tight-binding Hamiltonian (DFTB3-D3) to relax all the inter- and intra-molecular degrees of freedom of several thousands of generated crystal structures of five pharmaceutical-like molecules, saving a large amount of computational effort compared to earlier studies. The computational cost scales better with molecular size and flexibility than other CSP methods, suggesting that it could be extended to even larger and more flexible molecules. On average, this optimisation improved the average reproduction of the eight experimental crystal structures (RMSD15) and experimental conformers (RMSD1) by 4% and 23%, respectively. The intermolecular interactions were then further optimised using distributed multipoles, derived from the molecular wave-functions, to accurately describe the electrostatic components of the intermolecular energies. In all cases, the experimental crystal structures are close to the top of the lattice energy ranking. Phonon calculations on some of the lowest energy structures were also performed with DFTB3-D3 methods to calculate the vibrational component of the Helmholtz free energy, providing further insights into the solid-state behaviour of the target molecules. We conclude that DFTB3-D3 is a cost-effective method for optimising flexible molecules, bridging the gap between the approximate methods used in CSP searches for generating crystal structures and more accurate methods required in the final energy ranking.

Identifying pragmatic quasi-harmonic electronic structure approaches for modeling molecular crystal thermal expansion

Quasi-harmonic approaches provide an economical route to modeling the temperature dependence of molecular crystal structures and properties. Several studies have demonstrated good performance of these models, at least for rigid molecules, when using fragment-based approaches with correlated wavefunction techniques. Many others have found success employing dispersion-corrected density functional theory (DFT). Here, a hierarchy of models in which the energies, geometries, and phonons are computed either with correlated methods or DFT are examined to identify which combinations produce useful predictions for properties such as the molar volume, enthalpy, and entropy as a function of temperature. The results demonstrate that refining DFT geometries and phonons with single-point energies based on dispersion-corrected second-order Møller-Plesset perturbation theory can provide clear improvements in the molar volumes and enthalpies compared to those obtained from DFT alone. Predicted entropies, which are governed by vibrational contributions, benefit less clearly from the hybrid schemes. Using these hybrid techniques, the room-temperature thermochemistry of acetaminophen (paracetamol) is predicted to address the discrepancy between two experimental sublimation enthalpy measurements.

First-principles stability ranking of molecular crystal polymorphs with the DFT+MBD approach

The ability to accurately calculate the relative stabilities of numerous polymorphs of a given molecular crystal is crucial for the success of any molecular crystal structure prediction (CSP) approach. We have recently presented a hierarchical CSP procedure based on van-der-Waals-inclusive density functional theory [Hoja et al., 2018, arXiv:1803.07503], which yields excellent stability rankings for molecular crystals involving rigid molecules, salts, co-crystals, and highly polymorphic drug-like molecules. This approach includes many-body dispersion effects, exact exchange, as well as vibrational free energies. Here, we discuss in detail the impact of these effects on the obtained stability rankings. In addition, we assess the impact of the approximations used in our hierarchical procedure. We show that our procedure is generally robust to 1-2 kJ mol-1 for the systems in the latest CSP blind test but vibrational free energies for crystals involving flexible molecules would benefit from directly including many-body dispersion interactions. In addition, we also discuss the effect of temperature on the structure of molecular crystals and a simple but effective method for estimating anharmonic effects.

Thermal Expansion of Carbamazepine: Systematic Crystallographic Measurements Challenge Quantum Chemical Calculations

We report systematic temperature-dependent X-ray measurements on the most stable carbamazepine polymorph. This active pharmaceutical ingredient is used to demonstrate how the thermal expansion can probe certain intermolecular interactions resulting in anisotropic expansion behavior. We show that most structural features can be captured by electronic structure calculations at the quasi-harmonic approximation (QHA) provided a dispersion-corrected density functional based method is employed. The impact of thermal expansion on the phonon modes and hence free energy contributions is large enough to impact the relative stability of different polymorphs.

Application of computational methods to the design and characterisation of porous molecular materials

Composed from discrete units, porous molecular materials (PMMs) possess properties not observed for conventional, extended solids. Molecular simulations provide crucial understanding for the design and characterisation of these unique materials.