The potential of computed crystal energy landscapes to aid solid-form development

Dual space divergence in small molecule quasiracemates: benzoyl leucine and phenylalanine assemblies

Quasiracemic materials constructed with two points of structural difference were used to understand the role molecular shape plays in molecular assembly. Hot stage, crystallographic and occupied cavity space assessments provide insight into how imposed CH3/Cl and H/CF3 structural variations placed on benzoyl leucine and phenylalanine scaffolds result in a remarkably high occurrence of cocrystal formation.

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.

Predicting crystal form stability under real-world conditions

The physicochemical properties of molecular crystals, such as solubility, stability, compactability, melting behaviour and bioavailability, depend on their crystal form1. In silico crystal form selection has recently come much closer to realization because of the development of accurate and affordable free-energy calculations2-4. Here we redefine the state of the art, primarily by improving the accuracy of free-energy calculations, constructing a reliable experimental benchmark for solid-solid free-energy differences, quantifying statistical errors for the computed free energies and placing both hydrate crystal structures of different stoichiometries and anhydrate crystal structures on the same energy landscape, with defined error bars, as a function of temperature and relative humidity. The calculated free energies have standard errors of 1-2 kJ mol-1 for industrially relevant compounds, and the method to place crystal structures with different hydrate stoichiometries on the same energy landscape can be extended to other multi-component systems, including solvates. These contributions reduce the gap between the needs of the experimentalist and the capabilities of modern computational tools, transforming crystal structure prediction into a more reliable and actionable procedure that can be used in combination with experimental evidence to direct crystal form selection and establish control5.

A phase variety of fluorinated ionic liquids: Molecular conformational and crystal polymorph

Crystal polymorphs of fluorinated ionic liquids (fILs) were examined at low-temperature (LT) by Raman spectroscopy. The fILs were 1-alkyl-3-methylimidazolium perfluorobutanesulfonate, [C_nmim][PFBS] (n = 4, 6, and 8). The cations and anion possess conformational degrees of freedom. Various LT phases were derived from the conformational polymorphs of the cations and the anion. Conformational flexibility depended on alkyl chain length. The crystal polymorphs in the fILs were sensitive to molecular conformations and flexibility.

XDM-corrected hybrid DFT with numerical atomic orbitals predicts molecular crystal lattice energies with unprecedented accuracy

Molecular crystals are important for many applications, including energetic materials, organic semiconductors, and the development and commercialization of pharmaceuticals. The exchange-hole dipole moment (XDM) dispersion model has shown good performance in the calculation of relative and absolute lattice energies of molecular crystals, although it has traditionally been applied in combination with plane-wave/pseudopotential approaches. This has limited XDM to use with semilocal functional approximations, which suffer from delocalization error and poor quality conformational energies, and to systems with a few hundreds of atoms at most due to unfavorable scaling. In this work, we combine XDM with numerical atomic orbitals, which enable the efficient use of XDM-corrected hybrid functionals for molecular crystals. We test the new XDM-corrected functionals for their ability to predict the lattice energies of molecular crystals for the X23 set and 13 ice phases, the latter being a particularly stringent test. A composite approach using a XDM-corrected, 25% hybrid functional based on B86bPBE achieves a mean absolute error of 0.48 kcal mol^-1 per molecule for the X23 set and 0.19 kcal mol^-1 for the total lattice energies of the ice phases, compared to recent diffusion Monte-Carlo data. These results make the new XDM-corrected hybrids not only far more computationally efficient than previous XDM implementations, but also the most accurate density-functional methods for molecular crystal lattice energies to date.

Minimal Basis Set Hartree-Fock Corrected with Atom-Centered Potentials for Molecular Crystal Modeling and Crystal Structure Prediction

Crystal structure prediction (CSP), determining the experimentally observable structure of a molecular crystal from the molecular diagram, is an important challenge with technologically relevant applications in materials manufacturing and drug design. For the purpose of screening the randomly generated candidate crystal structures, CSP protocols require energy ranking methods that are fast and can accurately capture the small energy differences between molecular crystals. In addition, a good ranking method should also produce accurate equilibrium geometries, both intramolecular and intermolecular. In this article, we explore the combination of minimal-basis-set Hartree-Fock (HF) with atom-centered potentials (ACPs) as a method for modeling the structure and energetics of molecular crystals. The ACPs are developed for the H, C, N, and O atoms and fitted to a set of reference data at the B86bPBE-XDM level in order to mitigate basis-set incompleteness and missing correlation. In particular, ACPs are developed in combination with two methods: HF-D3/MINIs and HF-3c. The application of ACPs greatly improves the performance of HF-D3/MINIs for lattice energies, crystal energy differences, energy-volume and energy-strain relations, and crystal geometries. In the case of HF-3c, the improvement in the crystal energy differences is much smaller than in HF-D3/MINIs, but lattice energies and particularly crystal geometries are considerably better when ACPs are used. The resulting methods may be useful for CSP but also for quick calculation of molecular crystal lattice energies and geometries.

From Powders to Single Crystals: A Crystallographer's Toolbox for Small-Molecule Structure Determination

Although the crystal structures of small-molecule compounds are often determined from single-crystal X-ray diffraction (scXRD), recent advances in three-dimensional electron diffraction (3DED) and crystal structure prediction (CSP) methods promise to expand the structure elucidation toolbox available to the crystallographer. Herein, a comparative assessment of scXRD, 3DED, and CSP in combination with powder X-ray diffraction is carried out on two former drug candidate compounds and a multicomponent crystal of a key building block in the synthesis of gefapixant citrate.

Complete resonance assignment of a pharmaceutical drug at natural isotopic abundance from DNP-Enhanced solid-state NMR

Solid-state dynamic nuclear polarization enhanced magic angle spinning (DNP-MAS) NMR measurements coupled with density functional theory (DFT) calculations enable the full resonance assignment of a complex pharmaceutical drug molecule without the need for isotopic enrichment. DNP dramatically enhances the NMR signals, thereby making possible previously intractable two-dimensional correlation NMR spectra at natural abundance. Using inputs from DFT calculations, herein we describe a significant improvement to the structure elucidation process for complex organic molecules. Further, we demonstrate that a series of two-dimensional correlation experiments, including ¹⁵N-¹³C TEDOR, ¹³C-¹³C INADEQUATE/SARCOSY, ¹⁹F-¹³C HETCOR, and ¹H-¹³C HETCOR, can be obtained at natural isotopic abundance within reasonable experiment times, thus enabling a complete resonance assignment of sitagliptin, a pharmaceutical used for the treatment of type 2 diabetes.

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.

How many more polymorphs of ROY remain undiscovered

With 12 crystal forms, 5-methyl-2-[(2-nitrophenyl)amino]-3-thiophenecabonitrile (a.k.a. ROY) holds the current record for the largest number of fully characterized organic crystal polymorphs. Four of these polymorph structures have been reported since 2019, raising the question of how many more ROY polymorphs await future discovery. Employing crystal structure prediction and accurate energy rankings derived from conformational energy-corrected density functional theory, this study presents the first crystal energy landscape for ROY that agrees well with experiment. The lattice energies suggest that the seven most stable ROY polymorphs (and nine of the twelve lowest-energy forms) on the Z' = 1 landscape have already been discovered experimentally. Discovering any new polymorphs at ambient pressure will likely require specialized crystallization techniques capable of trapping metastable forms. At pressures above 10 GPa, however, a new crystal form is predicted to become enthalpically more stable than all known polymorphs, suggesting that further high-pressure experiments on ROY may be warranted. This work highlights the value of high-accuracy crystal structure prediction for solid-form screening and demonstrates how pragmatic conformational energy corrections can overcome the limitations of conventional density functionals for conformational polymorphs.

Along the road to Crystal Structure Prediction (CSP) of pharmaceutical-like molecules

Computational methods used for predicting crystal structures of organic compounds are mature enough to be routinely used with many rigid and semi-rigid organic molecules. The usefulness of Crystal Structure Prediction...

Solid-State Study of the Structure, Dynamics, and Thermal Processes of Safinamide Mesylate─A New Generation Drug for the Treatment of Neurodegenerative Diseases

Safinamide mesylate (SM), the pure active pharmaceutical ingredient (API) recently used in Parkinson disease treatment, recrystallized employing water–ethanol mixture of solvents (vol/vol 1:9) gives a different crystallographic form compared to SM in Xadago tablets. Pure SM crystallizes as a hemihydrate in the monoclinic system with the P2₁ space group. Its crystal and molecular structure were determined by means of cryo X-ray crystallography at 100 K. SM in the Xadago tablet exists in anhydrous form in the orthorhombic crystallographic system with the P2₁2₁2₁ space group. The water migration and thermal processes in the crystal lattice were monitored by solid-state NMR spectroscopy, differential scanning calorimetry, and thermogravimetric analysis. SM in Xadago in the high-humidity environment undergoes phase transformation to the P2₁ form which can be easily reversed just by heating up to 80 °C. For the commercial form of the API, there is also a reversible thermal transformation observed between Z′ = 1 ↔ Z′ = 3 crystallographic forms in the 0–20 °C temperature range. Analysis of molecular motion in the crystal lattice proves that the observed conformational polymorphism is forced by intramolecular dynamics. All above-mentioned processes were analyzed and described employing the NMR crystallography approach with the support of advanced theoretical calculations.

Expanding the Solid Form Landscape of Bipyridines

Two bipyridine isomers (2,2'- and 4,4'-), used as coformers and ligands in coordination chemistry, were subjected to solid form screening and crystal structure prediction. One anhydrate and a formic acid disolvate were crystallized for 2,2'-bipyridine, whereas multiple solid-state forms, anhydrate, dihydrate, and eight solvates with carboxylic acids, including a polymorphic acetic acid disolvate, were found for the 4,4'-isomer. Seven of the solvates are reported for the first time, and structural information is provided for six of the new solvates. All twelve solid-state forms were investigated comprehensively using experimental [thermal analysis, isothermal calorimetry, X-ray diffraction, gravimetric moisture (de)sorption, and IR spectroscopy] and computational approaches. Lattice and interaction energy calculations confirmed the thermodynamic driving force for disolvate formation, mediated by the absence of H-bond donor groups of the host molecules. The exposed location of the N atoms in 4,4'-bipyridine facilitates the accommodation of bigger carboxylic acids and leads to higher conformational flexibility compared to 2,2'-bipyridine. For the 4,4'-bipyridine anhydrate ↔ hydrate interconversion hardly any hysteresis and a fast transformation kinetics are observed, with the critical relative humidity being at 35% at room temperature. The computed anhydrate crystal energy landscapes have the 2,2'-bipyridine as the lowest energy structure and the 4,4'-bipyridine among the low-energy structures and suggest a different crystallization behavior of the two compounds.

Toward Physics-Based Solubility Computation for Pharmaceuticals to Rival Informatics

We demonstrate that physics-based calculations of intrinsic aqueous solubility can rival cheminformatics-based machine learning predictions. A proof-of-concept was developed for a physics-based approach via a sublimation thermodynamic cycle, building upon previous work that relied upon several thermodynamic approximations, notably the 2RT approximation, and limited conformational sampling. Here, we apply improvements to our sublimation free-energy model with the use of crystal phonon mode calculations to capture the contributions of the vibrational modes of the crystal. Including these improvements with lattice energies computed using the model-potential-based Ψ_mol method leads to accurate estimates of sublimation free energy. Combining these with hydration free energies obtained from either molecular dynamics free-energy perturbation simulations or density functional theory calculations, solubilities comparable to both experiment and informatics predictions are obtained. The application to coronene, succinic acid, and the pharmaceutical desloratadine shows how the methods must be adapted for the adoption of different conformations in different phases. The approach has the flexibility to extend to applications that cannot be covered by informatics methods.

Synergy of Solid-State NMR, Single-Crystal X-ray Diffraction, and Crystal Structure Prediction Methods: A Case Study of Teriflunomide (TFM)

In this work, for the first time, we present the X-ray diffraction crystal structure and spectral properties of a new, room-temperature polymorph of teriflunomide (TFM), CSD code 1969989. As revealed by DSC, the low-temperature TFM polymorph recently reported by Gunnam et al. undergoes a reversible thermal transition at -40 °C. This reversible process is related to a change in Z' value, from 2 to 1, as observed by variable-temperature ¹H-¹³C cross-polarization (CP) magic-angle spinning (MAS) solid-state NMR, while the crystallographic system is preserved (triclinic). Two-dimensional ¹³C-¹H and ¹H-¹H double-quantum MAS NMR spectra are consistent with the new room-temperature structure, including comparison with GIPAW (gauge-including projector augmented waves) calculated NMR chemical shifts. A crystal structure prediction procedure found both experimental teriflunomide polymorphs in the energetic global minimum region. Differences between the polymorphs are seen for the torsional angle describing the orientation of the phenyl ring relative to the planarity of the TFM molecule. In the low-temperature structure, there are two torsion angles of 4.5 and 31.9° for the two Z' = 2 molecules, while in the room-temperature structure, there is disorder that is modeled with ∼50% occupancy between torsion angles of -7.8 and 28.6°. These observations are consistent with a broad energy minimum as revealed by DFT calculations. PISEMA solid-state NMR experiments show a reduction in the C-H dipolar coupling in comparison to the static limit for the aromatic CH moieties of 75% and 51% at 20 and 40 °C, respectively, that is indicative of ring flips at the higher temperature. Our study shows the power of combining experiments, namely DSC, X-ray diffraction, and MAS NMR, with DFT calculations and CSP to probe and understand the solid-state landscape, and in particular the role of dynamics, for pharmaceutical molecules.

Modeling the α- and β-resorcinol phase boundary via combination of density functional theory and density functional tight-binding

The ability to predict not only what organic crystal structures might occur but also the thermodynamic conditions under which they are the most stable would be extremely useful for discovering and designing new organic materials. The present study takes a step in that direction by predicting the temperature- and pressure-dependent phase boundary between the α and β polymorphs of resorcinol using density functional theory (DFT) and the quasi-harmonic approximation. To circumvent the major computational bottleneck associated with computing a well-converged phonon density of states via the supercell approach, a recently developed approximation is employed, which combines a supercell phonon density of states from dispersion-corrected third-order density functional tight binding [DFTB3-D3(BJ)] with frequency corrections derived from a smaller B86bPBE-XDM functional DFT phonon calculation on the crystallographic unit cell. This mixed DFT/DFTB quasi-harmonic approach predicts the lattice constants and unit cell volumes to within 1%-2% at lower pressures. It predicts the thermodynamic phase boundary in almost perfect agreement with the experiment, although this excellent agreement does reflect fortuitous cancellation of errors between the enthalpy and entropy of transition.

The trimorphism of 3-hydroxybenzoic acid: an experimental and computational study

A computationally driven experimental search for polymorphs of 3-hydroxybenzoic acid confirmed the third form and the small energy differences between the polymorphs.

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'.

Free-energy landscape of polymer-crystal polymorphism

Polymorphism rationalizes how processing can control the final structure of a material. The rugged free-energy landscape and exceedingly slow kinetics in the solid state have so far hampered computational investigations. We report for the first time the free-energy landscape of a polymorphic crystalline polymer, syndiotactic polystyrene. Coarse-grained metadynamics simulations allow us to efficiently sample the landscape at large. The free-energy difference between the two main polymorphs, α and β, is further investigated by quantum-chemical calculations. The results of the two methods are in line with experimental observations: they predict β as the more stable polymorph under standard conditions. Critically, the free-energy landscape suggests how the α polymorph may lead to experimentally observed kinetic traps. The combination of multiscale modeling, enhanced sampling, and quantum-chemical calculations offers an appealing strategy to uncover complex free-energy landscapes with polymorphic behavior.

Virtual Coformer Screening by Crystal Structure Predictions: Crucial Role of Crystallinity in Pharmaceutical Cocrystallization

One of the most popular strategies of the optimization of drug properties in the pharmaceutical industry appears to be a solid form changing into a cocrystalline form. A number of virtual screening approaches have been previously developed to allow a selection of the most promising cocrystal formers (coformers) for an experimental follow-up. A significant drawback of those methods is related to the lack of accounting for the crystallinity contribution to cocrystal formation. To address this issue, we propose in this study two virtual coformer screening approaches based on a modern cloud-computing crystal structure prediction (CSP) technology at a dispersion-corrected density functional theory (DFT-D) level. The CSP-based methods were for the first time validated on challenging cases of indomethacin and paracetamol cocrystallization, for which the previously developed approaches provided poor predictions. The calculations demonstrated a dramatic improvement of the virtual coformer screening performance relative to the other methods. It is demonstrated that the crystallinity contribution to the formation of paracetamol and indomethacin cocrystals is a dominant one and, therefore, should not be ignored in the virtual screening calculations. Our results encourage a broad utilization of the proposed CSP-based technology in the pharmaceutical industry as the only virtual coformer screening method that directly accounts for the crystallinity contribution.

The Eight Hydrates of Strychnine Sulfate

Commercial samples of strychnine sulfate were used as the starting material in crystallization experiments accompanied by stability studies. Eight hydrate forms (HyA-HyG), including five novel hydrates, were verified. The crystal structures of HyA ("pentahydrate") and HyF ("hexahydrate") were determined from single-crystal X-ray diffraction data. HyF was identified as the most stable hydrate at high water activities at room temperature (RT), and HyA and HyC were also found to be stable at ambient conditions. Long-time storage experiments over nearly two decades confirm that these three hydrates are stable at ambient conditions (20-60% relative humidity). The other five hydrates, HyB ("dihydrate"), HyD, HyE, HyG, and HyH, are only observable at the low(est) relative humidity (RH) levels at RT. Some of these latter forms can only exist within a very narrow RH range and are therefore intermediate phases. By applying a range of complementary experimental techniques such as gravimetric moisture sorption analysis, thermal analysis, moisture controlled PXRD measurements, and variable temperature IR spectroscopy in combination with principal component analysis, it was possible to identify the distinct hydrate phases and elucidate their stability and dehydration pathways. The observed (de)hydration routes, HyA ↔ HyB, HyC ↔ HyD ↔ HyE, HyF ↔ HyG ↔ HyH and HyF → HyA ↔ HyB, depended on the initial hydrate form, particle size, and atmospheric conditions. In addition, a transformation from HyC/HyA to HyF occurs at high RH values at RT. The specific moisture and temperature conditions of none of the applied drying regimes yielded a crystalline water-free form, which highlights the essential role of water molecules for the formation and stability of the crystalline strychnine sulfate phases.

Calculation of Diamagnetic Susceptibility Tensors of Organic Crystals: From Coronene to Pharmaceutical Polymorphs

Understanding why crystallization in strong magnetic fields can lead to new polymorphs requires methods to calculate the diamagnetic response of organic molecular crystals. We develop the calculation of the macroscopic diamagnetic susceptibility tensor, χ^cryst, for organic molecular crystals using periodic density functional methods. The crystal magnetic susceptibility tensor, χ^cryst, for all experimentally known polymorphs, and its molecular counterpart, χ^mol, are calculated for flexible pharmaceuticals such as carbamazepine, flufenamic acid, and chalcones, and rigid molecules, such as benzene, pyridine, acridine, anthracene, and coronene, whose molecular magnetic properties have been traditionally studied. A tensor addition method is developed to approximate the crystal diamagnetic susceptibility tensor, χ^cryst, from the molecular one, χ^mol, giving good agreement with those calculated directly using the more costly periodic density functional method for χ^cryst. The response of pharmaceutical molecules and crystals to magnetic fields, as embodied by χ^cryst, is largely determined by the packing in the crystal, as well as the molecular conformation. The anisotropy of χ^cryst can vary considerably between polymorphs though the isotropic terms are fairly constant. The implications for developing a computational method for predicting whether crystallization in a magnetic field could produce a novel or different polymorph are discussed.

Charge mobility calculation of organic semiconductors without use of experimental single-crystal data

Prediction of material properties of newly designed molecules is a long-term goal in organic electronics. In general, it is a difficult problem, because the material properties are dominated by the unknown packing structure. We present a practical method to obtain charge transport properties of organic single crystals, without use of experimental single-crystal data. As a demonstration, we employ the promising molecule C10-DNBDT. We succeeded in quantitative evaluation of charge mobility of the single crystal using our quantum wave-packet dynamical simulation method. Here, the single-crystal data is computationally obtained by searching possible packing structures from structural formula of the molecule. We increase accuracy in identifying the actual crystal structure from suggested ones by using not only crystal energy but also similarity between calculated and experimental powder X-ray diffraction patterns. The proposed methodology can be a theoretical design technique for efficiently developing new high-performance organic semiconductors, since it can estimate the charge transport properties at early stage in the process of material development.

Fast energy minimization of the CCDC drug-subset structures by molecule-in-cluster computations allows independent structure validation and model completion

Optimizing structures with computations on clusters of molecules permits generation of structure-specific restraints for refinement and structure validation.

Experimental and computational approaches to rationalise multicomponent supramolecular assemblies: dapsone monosolvates

The monosolvate crystal energy landscapes of dapsone (DDS) including the solvents carbon tetrachloride, acetone, cyclohexanone, dimethyl formamide, tetrahydrofuran, methyl ethyl ketone, 1,2-dichloroethane, 1,4-dioxane, dichloromethane and chloroform were established using experimental and computational approaches. To rationalise and understand solvate formation, solvate stability and desolvation reactions a careful control of the experimental crystallisation and storage conditions, a range of thermoanalytical methods and crystal structure prediction were required. Six of the eight DDS monosolvates are reported and characterised for the first time. Structural similarity and diversity of the at ambient conditions unstable monosolvates were apparent from the computed crystal energy landscapes, which had the experimental packings as lowest energy structures. The computed structures were used as input for Rietveld refinements and isostructurality of four of the monosolvates was confirmed. Packing comparisons of the solvate structures and molecular properties of the solvent molecules indicated that both size/shape of the solvent molecule and the possible DDSsolvent interactions are the important factors for DDS solvate formation. Through the combination of experiment and theory solvate stability and structural features have been rationalised.

Charge Distributions of Nitro Groups Within Organic Explosive Crystals: Effects on Sensitivity and Modeling

The charge distribution of NO₂ groups within the crystalline polymorphs of energetic materials strongly affects their explosive properties. We use the recently introduced basis-space iterated stockholder atom partitioning of high-quality charge distributions to examine the approximations that can be made in modeling polymorphs and their physical properties, using 1,3,5-trinitroperhydro-1,3,5-triazine, trinitrotoluene, 1-3-5-trinitrobenzene, and hexanitrobenzene as exemplars. The NO₂ charge distribution is strongly affected by the neighboring atoms, the rest of the molecules, and also significantly by the NO₂ torsion angle within the possible variations found in observed crystal structures. Thus, the proposed correlations between the molecular electrostatic properties, such as trigger-bond potential or maxima in the electrostatic potential, and impact sensitivity will be affected by the changes in conformation that occur on crystallization. We establish the relationship between the NO₂ torsion angle and the likelihood of occurrence in observed crystal structures, the conformational energy, and the charge and dipole magnitude on each atom, and how this varies with the neighboring groups. We examine the effect of analytically rotating the atomic multipole moments to model changes in torsion angle and establish that this is a viable approach for crystal structures but is not accurate enough to model the relative lattice energies. This establishes the basis of transferability of the NO₂ charge distribution for realistic nonempirical model intermolecular potentials for simulating energetic materials.

Dapsone Form V: A Late Appearing Thermodynamic Polymorph of a Pharmaceutical

Five anhydrate polymorphs (forms I-V) and the isomorphic dehydrate (Hy_dehy) of dapsone (4,4'-diaminodiphenyl sulfone or DDS) were prepared and characterized in an interdisciplinary experimental and computational study, elucidating the kinetic and thermodynamic stabilities, solid form interrelationships, and structural features of the known forms I-IV, the novel polymorph form V, and Hy_dehy. Calorimetric measurements, solubility experiments, and lattice energy calculations revealed that form V is the thermodynamically stable polymorph from absolute zero to at least 90 °C. At higher temperatures, form II, and then form I, becomes the most stable DDS solid form. The computed 0 K stability order (lattice energy calculations) was confirmed with calorimetric measurements as follows, V (most stable) > III > Hy_dehy > II > I > IV (least stable). The discovery of form V was complicated by the fact that the metastable but kinetically stabilized form III shows a higher nucleation and growth rate. By combining laboratory powder X-ray diffraction data and ab initio calculations, the crystal structure of form V ( P2₁/ c, Z' = 4) was solved, with a high energy DDS conformation allowing a denser packing and more stable intermolecular interactions, rationalizing the formation of a high Z' structure. The structures of the forms I and IV, only observed from the melt and showing distinct packing features compared to the forms II, III, and V, were derived from the computed crystal energy landscapes. Dehydration modeling of the DDS hydrate led to the Hy_dehy structure. This study expands our understanding about the complex crystallization behavior of pharmaceuticals and highlights the big challenge in solid form screening, especially that there is no clear end point.

The use of variable temperature 13 C solid-state MAS NMR and GIPAW DFT calculations to explore the dynamics of diethylcarbamazine citrate

Experimental ¹³ C solid-state magic-angle spinning (MAS) Nuclear Magnetic Resonance (NMR) as well as Density-Functional Theory (DFT) gauge-including projector augmented wave (GIPAW) calculations were used to probe disorder and local mobility in diethylcarbamazine citrate, (DEC)⁺ (citrate)^- . This compound has been used as the first option drug for the treatment of filariasis, a disease endemic in tropical countries and caused by adult worms of Wuchereria bancrofti, which is transmitted by mosquitoes. We firstly present 2D ¹³ C─¹ H dipolar-coupling-mediated heteronuclear correlation spectra recorded at moderate spinning frequency, to explore the intermolecular interaction between DEC and citrate molecules. Secondly, we investigate the dynamic behavior of (DEC)⁺ (citrate)^- by varying the temperature and correlating the experimental MAS NMR results with DFT GIPAW calculations that consider two (DEC)⁺ conformers (in a 70:30 ratio) for crystal structures determined at 293 and 235 K. Solid-state NMR provides insights on slow exchange dynamics revealing conformational changes involving particularly the DEC ethyl groups.

A generally applicable atomic-charge dependent London dispersion correction

The so-called D4 model is presented for the accurate computation of London dispersion interactions in density functional theory approximations (DFT-D4) and generally for atomistic modeling methods. In this successor to the DFT-D3 model, the atomic coordination-dependent dipole polarizabilities are scaled based on atomic partial charges which can be taken from various sources. For this purpose, a new charge-dependent parameter-economic scaling function is designed. Classical charges are obtained from an atomic electronegativity equilibration procedure for which efficient analytical derivatives with respect to nuclear positions are developed. A numerical Casimir-Polder integration of the atom-in-molecule dynamic polarizabilities then yields charge- and geometry-dependent dipole-dipole dispersion coefficients. Similar to the D3 model, the dynamic polarizabilities are precomputed by time-dependent DFT and all elements up to radon (Z = 86) are covered. The two-body dispersion energy expression has the usual sum-over-atom-pairs form and includes dipole-dipole as well as dipole-quadrupole interactions. For a benchmark set of 1225 molecular dipole-dipole dispersion coefficients, the D4 model achieves an unprecedented accuracy with a mean relative deviation of 3.8% compared to 4.7% for D3. In addition to the two-body part, three-body effects are described by an Axilrod-Teller-Muto term. A common many-body dispersion expansion was extensively tested, and an energy correction based on D4 polarizabilities is found to be advantageous for larger systems. Becke-Johnson-type damping parameters for DFT-D4 are determined for more than 60 common density functionals. For various standard energy benchmark sets, DFT-D4 slightly but consistently outperforms DFT-D3. Especially for metal containing systems, the introduced charge dependence of the dispersion coefficients improves thermochemical properties. We suggest (DFT-)D4 as a physically improved and more sophisticated dispersion model in place of DFT-D3 for DFT calculations as well as other low-cost approaches like semi-empirical models.

Crystal-energy landscapes of active pharmaceutical ingredients using composite approaches

Composite methods employing dispersion-corrected DFT consistently identify experimentally isolated polymorphs as the lowest-energy crystal structures of common APIs.

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.

Energy-Structure-Function Maps: Cartography for Materials Discovery

Some of the most successful approaches to structural design in materials chemistry have exploited strong directional bonds, whose geometric reliability lends predictability to solid-state assembly. For example, metal-organic frameworks are an important design platform in materials chemistry. By contrast, the structure of molecular crystals is defined by a balance of weaker intermolecular forces, and small changes to the molecular building blocks can lead to large changes in crystal packing. Hence, empirical rules are inherently less reliable for engineering the structures of molecular solids. Energy-structure-function (ESF) maps are a new approach for the discovery of functional organic crystals. These maps fuse crystal-structure prediction with the computation of physical properties to allow researchers to choose the most promising molecule for a given application, prior to its synthesis. ESF maps were used recently to discover a highly porous molecular crystal that has a high methane deliverable capacity and the lowest density molecular crystal reported to date (r = 0.41 g cm^-3 , SA_BET = 3425 m² g^-1 ). Progress in this field is reviewed, with emphasis on the future opportunities and challenges for a design strategy based on computed ESF maps.

Simplified DFT methods for consistent structures and energies of large systems

Kohn-Sham density functional theory (DFT) is routinely used for the fast electronic structure computation of large systems and will most likely continue to be the method of choice for the generation of reliable geometries in the foreseeable future. Here, we present a hierarchy of simplified DFT methods designed for consistent structures and non-covalent interactions of large systems with particular focus on molecular crystals. The covered methods are a minimal basis set Hartree-Fock (HF-3c), a small basis set screened exchange hybrid functional (HSE-3c), and a generalized gradient approximated functional evaluated in a medium-sized basis set (B97-3c), all augmented with semi-classical correction potentials. We give an overview on the methods design, a comprehensive evaluation on established benchmark sets for geometries and lattice energies of molecular crystals, and highlight some realistic applications on large organic crystals with several hundreds of atoms in the primitive unit cell.

Supramolecular Organization of Nonstoichiometric Drug Hydrates: Dapsone

The observed moisture- and temperature dependent transformations of the dapsone (4,4'-diaminodiphenyl sulfone, DDS) 0. 33-hydrate were correlated to its structure and the number and strength of the water-DDS intermolecular interactions. A combination of characterization techniques was used, including thermal analysis (hot-stage microscopy, differential scanning calorimetry and thermogravimetric analysis), gravimetric moisture sorption/desorption studies and variable humidity powder X-ray diffraction, along with computational modeling (crystal structure prediction and pair-wise intermolecular energy calculations). Depending on the relative humidity the hydrate contains between 0 and 0.33 molecules of water per molecule DDS. The crystal structure is retained upon dehydration indicating that DDS hydrate shows a non-stoichiometric (de)hydration behavior. Unexpectedly, the water molecules are not located in structural channels but at isolated-sites of the host framework, which is counterintuitively for a hydrate with non-stoichiometric behavior. The water-DDS interactions were estimated to be weaker than water-host interactions that are commonly observed in stoichiometric hydrates and the lattice energies of the isomorphic dehydration product (hydrate structure without water molecules) and (form III) differ only by ~1 kJ mol-1. The computational generation of hypothetical monohydrates confirms that the hydrate with the unusual DDS:water ratio of 3:1 is more stable than a feasible monohydrate structure. Overall, this study highlights that a deeper understanding of the formation of hydrates with non-stoichiometric behavior requires a multidisciplinary approach including suitable experimental and computational methods providing a firm basis for the development and manufacturing of high quality drug products.

Fast and accurate quantum Monte Carlo for molecular crystals

Computational approaches based on the fundamental laws of quantum mechanics are now integral to almost all materials design initiatives in academia and industry. If computational materials science is genuinely going to deliver on its promises, then an electronic structure method with consistently high accuracy is urgently needed. We show that, thanks to recent algorithmic advances and the strategy developed in our manuscript, quantum Monte Carlo yields extremely accurate predictions for the lattice energies of materials at a surprisingly modest computational cost. It is thus no longer a technique that requires a world-leading computational facility to obtain meaningful results. While we focus on molecular crystals, the significance of our findings extends to all classes of materials.

Computer simulation plays a central role in modern-day materials science. The utility of a given computational approach depends largely on the balance it provides between accuracy and computational cost. Molecular crystals are a class of materials of great technological importance which are challenging for even the most sophisticated ab initio electronic structure theories to accurately describe. This is partly because they are held together by a balance of weak intermolecular forces but also because the primitive cells of molecular crystals are often substantially larger than those of atomic solids. Here, we demonstrate that diffusion quantum Monte Carlo (DMC) delivers subchemical accuracy for a diverse set of molecular crystals at a surprisingly moderate computational cost. As such, we anticipate that DMC can play an important role in understanding and predicting the properties of a large number of molecular crystals, including those built from relatively large molecules which are far beyond reach of other high-accuracy methods.

Structure and Stability of Molecular Crystals with Many-Body Dispersion-Inclusive Density Functional Tight Binding

Accurate prediction of structure and stability of molecular crystals is crucial in materials science and requires reliable modeling of long-range dispersion interactions. Semiempirical electronic structure methods are computationally more efficient than their ab initio counterparts, allowing structure sampling with significant speedups. We combine the Tkatchenko-Scheffler van der Waals method (TS) and the many-body dispersion method (MBD) with third-order density functional tight-binding (DFTB3) via a charge population-based method. We find an overall good performance for the X23 benchmark database of molecular crystals, despite an underestimation of crystal volume that can be traced to the DFTB parametrization. We achieve accurate lattice energy predictions with DFT+MBD energetics on top of vdW-inclusive DFTB3 structures, resulting in a speedup of up to 3000 times compared with a full DFT treatment. This suggests that vdW-inclusive DFTB3 can serve as a viable structural prescreening tool in crystal structure prediction.

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.

A variable-temperature X-ray diffraction and theoretical study of conformational polymorphism in a complex organic molecule (DTC)

Two conformational crystal polymorphs of 3-diethylamino-4-(4-methoxyphenyl)-1,1-dioxo-4H-1λ⁶,2-thiazete-4-carbonitrile (DTC) have been analyzed in the 100 K-room temperature range by single crystal X-ray diffraction and high quality DFT calculations. DTC has strongly polar as well as aliphatic substituents but no hydrogen bonding groups, and thus qualifies as a test molecule for the relative importance of electrostatic vs. dispersion–repulsion terms. The two polymorphs have the same P2₁/n space group and differ by a flipping of the –OCH₃ group, the two conformations being almost equi-energetic and separated by a low barrier. The system is monotropic in the observed temperature range with nearly identical thermal expansion coefficients and energy–temperature slopes, one phase consistently predicted to be more stable in agreement with the relative ease of appearance. Energy decompositions show that the electrostatic term is dominant and stabilizes with decreasing temperature. Dispersion and repulsion show the expected behavior, the former becoming more stabilizing at lower temperature in contrast with increasing repulsion at higher density. Absolute values and trends are very similar in the two phases, explaining the small total energy difference. Geometrical analyses of intermolecular contacts using fingerprint plots, as well as the study of molecular dipole moments as a function of T in the framework of the Quantum Theory of Atoms in Molecules, reveal more details of phase stability.

A small conformational change in the asymmetric unit has a significant effect on how non-covalent interactions determine (i) the crystal packing and (ii) the effect of T on the relative balance of electrostatics and dispersion–repulsions.

Is zeroth order crystal structure prediction (CSP_0) coming to maturity? What should we aim for in an ideal crystal structure prediction code?

Crystal structure prediction based on searching for the global minimum in the lattice energy (CSP_0) is growing in use for guiding the discovery of new materials, for example, new functional materials, new phases of interest to planetary scientists and new polymorphs relevant to pharmaceutical development. This Faraday Discussion can assess the progress of CSP_0 over the range of types of materials to which CSP is currently and could be applied, which depends on our ability to model the variety of interatomic forces in crystals. The basic hypothesis, that the outcome of crystallisation is determined by thermodynamics, needs examining by considering methods of modelling relative thermodynamic stability not only as a function of pressure and temperature, but also of size, solvent and the presence of heterogeneous templates or impurities (CSP_thd). Given that many important materials persist, and indeed may be formed, when they are not the most thermodynamically stable structure, we need to define what would be required of an ideal CSP code (CSP_aim).

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.

Commentary: Why Pharmaceutical Scientists in Early Drug Discovery Are Critical for Influencing the Design and Selection of Optimal Drug Candidates

This commentary reflects the collective view of pharmaceutical scientists from four different organizations with extensive experience in the field of drug discovery support. Herein, engaging discussion is presented on the current and future approaches for the selection of the most optimal and developable drug candidates. Over the past two decades, developability assessment programs have been implemented with the intention of improving physicochemical and metabolic properties. However, the complexity of both new drug targets and non-traditional drug candidates provides continuing challenges for developing formulations for optimal drug delivery. The need for more enabled technologies to deliver drug candidates has necessitated an even more active role for pharmaceutical scientists to influence many key molecular parameters during compound optimization and selection. This enhanced role begins at the early in vitro screening stages, where key learnings regarding the interplay of molecular structure and pharmaceutical property relationships can be derived. Performance of the drug candidates in formulations intended to support key in vivo studies provides important information on chemotype-formulation compatibility relationships. Structure modifications to support the selection of the solid form are also important to consider, and predictive in silico models are being rapidly developed in this area. Ultimately, the role of pharmaceutical scientists in drug discovery now extends beyond rapid solubility screening, early form assessment, and data delivery. This multidisciplinary role has evolved to include the practice of proactively taking part in the molecular design to better align solid form and formulation requirements to enhance developability potential.

From serendipity to supramolecular design: assessing the utility of computed crystal form landscapes in inferring the risks of crystal hydration in carboxylic acids

Calculated structural descriptors for predicted anhydrate polymorphs are used to assess the risks of crystal hydration in carboxylic acids.

Towards the systematic crystallisation of molecular ionic cocrystals: insights from computed crystal form landscapes

The underlying molecular and crystal properties affecting the crystallisation of organic molecular ionic cocrystals (ICCs) are investigated.