Report on the sixth blind test of organic crystal structure prediction methods

The application of tailor-made force fields and molecular dynamics for NMR crystallography: a case study of free base cocaine

Motional averaging has been proven to be significant in predicting the chemical shifts in ab initio solid-state NMR calculations, and the applicability of motional averaging with molecular dynamics has been shown to depend on the accuracy of the molecular mechanical force field. The performance of a fully automatically generated tailor-made force field (TMFF) for the dynamic aspects of NMR crystallography is evaluated and compared with existing benchmarks, including static dispersion-corrected density functional theory calculations and the COMPASS force field. The crystal structure of free base cocaine is used as an example. The results reveal that, even though the TMFF outperforms the COMPASS force field for representing the energies and conformations of predicted structures, it does not give significant improvement in the accuracy of NMR calculations. Further studies should direct more attention to anisotropic chemical shifts and development of the method of solid-state NMR calculations.

Ab initio crystal structure prediction of magnesium (poly)sulfides and calculation of their NMR parameters

Ab initio prediction of sensible crystal structures can be regarded as a crucial task in the quickly-developing methodology of NMR crystallography. In this contribution, an evolutionary algorithm was used for the prediction of magnesium (poly)sulfide crystal structures with various compositions. The employed approach successfully identified all three experimentally detected forms of MgS, i.e. the stable rocksalt form and the metastable wurtzite and zincblende forms. Among magnesium polysulfides with a higher content of sulfur, the most probable structure with the lowest formation energy was found to be MgS2, exhibiting a modified rocksalt structure, in which S2− anions were replaced by S2 2− dianions. Magnesium polysulfides with even larger fractions of sulfur were not predicted to be stable. For the lowest-energy structures, 25Mg quadrupolar coupling constants and chemical shift parameters were calculated using the density functional theory approach. The calculated NMR parameters could be well rationalized by the symmetries of the local magnesium environments, by the coordination of magnesium cations and by the nature of the surrounding anions. In the future, these parameters could serve as a reference for the experimentally determined 25Mg NMR parameters of magnesium sulfide species.

Dynamic quantum crystallography: lattice-dynamical models refined against diffraction data. II. Applications to L-alanine, naphthalene and xylitol

In the first paper of this series [Hoser & Madsen (2016). Acta Cryst. A72, 206-214], a new approach was introduced which enables the refinement of frequencies of normal modes obtained from ab initio periodic computations against single-crystal diffraction data. In this contribution, the performance of this approach is tested by refinement against data in the temperature range from 23 to 205 K on the molecular crystals of L-alanine, naphthalene and xylitol. The models, which are lattice-dynamical models derived at the Γ point of the Brillouin zone, are able to describe the atomic vibrations of L-alanine and naphthalene to a level where the residual densities are similar to those obtained from the independent atom model. For the more flexible molecule xylitol, larger deviations are found. Hydrogen ADPs (anisotropic displacement parameters) derived from the models are in similar or better agreement with neutron diffraction results than ADPs obtained by other procedures. The heat capacity calculated after normal mode refinement for naphthalene is in reasonable agreement with the heat capacity obtained from calorimetric measurements (to less than 1 cal mol^-1 K^-1 below 300 K), with deviations at higher temperatures indicating anharmonicity. Standard uncertainties and correlation of the refined parameters have been derived based on a Monte Carlo procedure. The uncertainties are quite small and probably underestimated.

Exchange-Hole Dipole Dispersion Model for Accurate Energy Ranking in Molecular Crystal Structure Prediction

Accurate energy ranking is a key facet to the problem of first-principles crystal-structure prediction (CSP) of molecular crystals. This work presents a systematic assessment of B86bPBE-XDM, a semilocal density functional combined with the exchange-hole dipole moment (XDM) dispersion model, for energy ranking using 14 compounds from the first five CSP blind tests. Specifically, the set of crystals studied comprises 11 rigid, planar compounds and 3 co-crystals. The experimental structure was correctly identified as the lowest in lattice energy for 12 of the 14 total crystals. One of the exceptions is 4-hydroxythiophene-2-carbonitrile, for which the experimental structure was correctly identified once a quasi-harmonic estimate of the vibrational free-energy contribution was included, evidencing the occasional importance of thermal corrections for accurate energy ranking. The other exception is an organic salt, where charge-transfer error (also called delocalization error) is expected to cause the base density functional to be unreliable. Provided the choice of base density functional is appropriate and an estimate of temperature effects is used, XDM-corrected density-functional theory is highly reliable for the energetic ranking of competing crystal structures.

Architecture of packing in molecular crystals

Molecular coordination shells have been used to classify crystal structures into similarity groups.

Accurate and efficient representation of intramolecular energy in ab initio generation of crystal structures. I. Adaptive local approximate models

The global search stage of crystal structure prediction (CSP) methods requires a fine balance between accuracy and computational cost, particularly for the study of large flexible molecules. A major improvement in the accuracy and cost of the intramolecular energy function used in the CrystalPredictor II [Habgood et al. (2015). J. Chem. Theory Comput. 11, 1957-1969] program is presented, where the most efficient use of computational effort is ensured via the use of adaptive local approximate model (LAM) placement. The entire search space of the relevant molecule's conformations is initially evaluated using a coarse, low accuracy grid. Additional LAM points are then placed at appropriate points determined via an automated process, aiming to minimize the computational effort expended in high-energy regions whilst maximizing the accuracy in low-energy regions. As the size, complexity and flexibility of molecules increase, the reduction in computational cost becomes marked. This improvement is illustrated with energy calculations for benzoic acid and the ROY molecule, and a CSP study of molecule (XXVI) from the sixth blind test [Reilly et al. (2016). Acta Cryst. B72, 439-459], which is challenging due to its size and flexibility. Its known experimental form is successfully predicted as the global minimum. The computational cost of the study is tractable without the need to make unphysical simplifying assumptions.

Crystal structures of p-substituted derivatives of 2,6-di-methyl-bromo-benzene with ½ ≤ Z' ≤ 4

The crystal structures of four bromo-arenes based on 2,6-di-methyl-bromo-benzene are reported, which are differentiated according the functional group X placed para to the Br atom: X = CN (4-bromo-3,5-di-methyl-benzo-nitrile, C9H8BrN), (1), X = NO2 (2-bromo-1,3-dimethyl-5-nitro-benzene, C8H8BrNO2), (2), X = NH2 (4-bromo-3,5-di-methyl-aniline, C8H10BrN), (3) and X = OH (4-bromo-3,5-di-methyl-phenol, C8H9BrO), (4). The content of the asymmetric unit is different in each crystal, Z' = ½ (X = CN), Z' = 1 (X = NO2), Z' = 2 (X = NH2), and Z' = 4 (X = OH), and is related to the mol-ecular symmetry and the propensity of X to be involved in hydrogen bonding. In none of the studied compounds does the crystal structure feature other non-covalent inter-actions, such as π-π, C-H⋯π or C-Br⋯Br contacts.

Modelling temperature-dependent properties of polymorphic organic molecular crystals

We present a large-scale study of the temperature-dependence of structures, free energy differences and properties of polymorphic molecular organic crystals. Lattice-vibrational Gibbs free energy differences between 475 pairs of polymorphs of organic molecular crystals have been calculated at 0 K and at their respective melting points, using a highly accurate anisotropic multipole-based force field and including thermal expansion through the use of a (negative) thermal pressure. Re-ranking of the relative thermodynamic stability of the polymorphs in each pair indicates the possibility of an enantiotropic phase transition between the crystal structures, which occurs in 21% of the studied systems. While vibrational contributions to free energies can have a significant effect on thermodynamic stability, the impact of thermal expansion on polymorph free energy differences is generally very small. We also calculate thermal expansion coefficients for the 864 crystal structures and investigate the temperature-dependence of mechanical properties, and pairwise differences in these properties between polymorphs.

Predicting Molecular Crystal Properties from First Principles: Finite-Temperature Thermochemistry to NMR Crystallography

Molecular crystals occur widely in pharmaceuticals, foods, explosives, organic semiconductors, and many other applications. Thanks to substantial progress in electronic structure modeling of molecular crystals, attention is now shifting from basic crystal structure prediction and lattice energy modeling toward the accurate prediction of experimentally observable properties at finite temperatures and pressures. This Account discusses how fragment-based electronic structure methods can be used to model a variety of experimentally relevant molecular crystal properties. First, it describes the coupling of fragment electronic structure models with quasi-harmonic techniques for modeling the thermal expansion of molecular crystals, and what effects this expansion has on thermochemical and mechanical properties. Excellent agreement with experiment is demonstrated for the molar volume, sublimation enthalpy, entropy, and free energy, and the bulk modulus of phase I carbon dioxide when large basis second-order Møller-Plesset perturbation theory (MP2) or coupled cluster theories (CCSD(T)) are used. In addition, physical insight is offered into how neglect of thermal expansion affects these properties. Zero-point vibrational motion leads to an appreciable expansion in the molar volume; in carbon dioxide, it accounts for around 30% of the overall volume expansion between the electronic structure energy minimum and the molar volume at the sublimation point. In addition, because thermal expansion typically weakens the intermolecular interactions, neglecting thermal expansion artificially stabilizes the solid and causes the sublimation enthalpy to be too large at higher temperatures. Thermal expansion also frequently weakens the lower-frequency lattice phonon modes; neglecting thermal expansion causes the entropy of sublimation to be overestimated. Interestingly, the sublimation free energy is less significantly affected by neglecting thermal expansion because the systematic errors in the enthalpy and entropy cancel somewhat. Second, because solid state nuclear magnetic resonance (NMR) plays an increasingly important role in molecular crystal studies, this Account discusses how fragment methods can be used to achieve higher-accuracy chemical shifts in molecular crystals. Whereas widely used plane wave density functional theory models are largely restricted to generalized gradient approximation (GGA) functionals like PBE in practice, fragment methods allow the routine use of hybrid density functionals with only modest increases in computational cost. In extensive molecular crystal benchmarks, hybrid functionals like PBE0 predict chemical shifts with 20-30% higher accuracy than GGAs, particularly for ¹H, ¹³C, and ¹⁵N nuclei. Due to their higher sensitivity to polarization effects, ¹⁷O chemical shifts prove slightly harder to predict with fragment methods. Nevertheless, the fragment model results are still competitive with those from GIPAW. The improved accuracy achievable with fragment approaches and hybrid density functionals increases discrimination between different potential assignments of individual shifts or crystal structures, which is critical in NMR crystallography applications. This higher accuracy and greater discrimination are highlighted in application to the solid state NMR of different acetaminophen and testosterone crystal forms.

Why do Hydrates (Solvates) Form in Small Neutral Organic Molecules? Exploring the Crystal Form Landscapes of the Alkaloids Brucine and Strychnine

Computational methods were used to generate and explore the crystal structure landscapes of the two alkaloids strychnine and brucine. The computed structures were analyzed and rationalized by correlating the modelling results to a rich pool of available experimental data. Despite their structural similarity, the two compounds show marked differences in the formation of solid forms. For strychnine only one anhydrous form is reported in the literature and two new solvates from 1,4-dioxane were detected in the course of this work. In contrast, 22 solid forms are so far known to exist for brucine, comprising two anhydrates, four hydrates (HyA - HyC and a 5.25-hydrate), twelve solvates (alcohols and acetone) and four heterosolvates (mixed solvates with water and alcohols). For strychnine it is hard to produce any solid form other than the stable anhydrate while the formation of specific solid state forms of brucine is governed by a complex interplay between temperature and relative humidity/water activity and it is rather a challenging to avoid hydrate formation. Differences in crystal packing and the high tendency for brucine to form hydrates are not intuitive from the molecular structure alone, as both molecules have hydrogen bond acceptor groups but lack hydrogen bond donor groups. Only the evaluation of the crystal energy landscapes, in particular the close-packed crystal structures and high-energy open frameworks containing voids of molecular (water) dimensions, allowed us to unravel the diverse solid state behavior of the two alkaloids at a molecular level. In this study we demonstrate that expanding the analysis of anhydrate crystal energy landscapes to higher energy structures and calculating the solvent-accessible volume can be used to estimate non-stoichiometric or channel hydrate (solvate) formation, without explicitly computing the hydrate/solvate crystal energy landscapes.

Solvent inclusion in the crystal structure of bis-[(adamantan-1-yl)methanaminium chloride] 1,4-dioxane hemisolvate monohydrate explained using the computed crystal energy landscape

Repeated attempts to crystallize 1-adamantane-methyl-amine hydro-chloride as an anhydrate failed but the salt was successfully crystallized as a solvate (2C11H20N+·2Cl-·0.5C4H8O2·H2O), with water and 1,4-dioxane playing a structural role in the crystal and engaging in hydrogen-bonding inter-actions with the cation and anion. Computational crystal-structure prediction was used to rationalize the solvent-inclusion behaviour of this salt by computing the solvent-accessible voids in the predicted low-energy structures for the anhydrate: the global lattice-energy minimum structure, which has the same packing of the ions as the solvate, has solvent-accessible voids that account for 3.71% of the total unit-cell volume and is 6 kJ mol-1 more stable than the next most stable predicted structure.

Crystal structure prediction: are we there yet?

This contribution comments on the advances of the latest Crystal Structure Prediction blind test and the challenges still lying ahead.

Effect of packing motifs on the energy ranking and electronic properties of putative crystal structures of tricyano-1,4-dithiino[c]-isothiazole

We present an analysis of putative structures of tricyano-1,4-dithiino[c]-isothiazole (TCS3), generated within the sixth crystal structure prediction blind test. Typical packing motifs are identified and characterized in terms of distinct patterns of close contacts and regions of electrostatic and dispersion interactions. We find that different dispersion-inclusive density functional theory (DFT) methods systematically favor specific packing motifs, which may affect the outcome of crystal structure prediction efforts. The effect of crystal packing on the electronic and optical properties of TCS3 is investigated using many-body perturbation theory within the GW approximation and the Bethe-Salpeter equation (BSE). We find that a structure with Pna21 symmetry and a bilayer packing motif exhibits intermolecular bonding patterns reminiscent of π-π stacking and has markedly different electronic and optical properties than the experimentally observed P21/n structure with a cyclic dimer motif, including a narrower band gap, enhanced band dispersion and broader optical absorption. The Pna21 bilayer structure is close in energy to the observed structure and may be feasible to grow.

Crystal structure prediction of rigid molecules

A non-polarizable force field based on atomic multipoles fit to reproduce experimental crystal properties andab initiogas-phase dimers is described. The Ewald method is used to calculate both long-range electrostatic and 1/r6dispersion energies of crystals. The dispersion energy of a crystal calculated by a cutoff method is shown to converge slowly to the exact Ewald result. A method for constraining space-group symmetry during unit-cell optimization is derived. Results for locally optimizing 4427 unit cells including volume, cell parameters, unit-cell r.m.s.d. and CPU timings are given for both flexible and rigid molecule optimization. An algorithm for randomly generating rigid molecule crystals is described. Using the correct experimentally determined space group, the average and maximum number of random crystals needed to find the correct experimental structure is given for 2440 rigid single component crystals. The force field energy rank of the correct experimental structure is presented for the same set of 2440 rigid single component crystals assuming the correct space group. A complete crystal prediction is performed for two rigid molecules by searching over the 32 most probable space groups.

Exploring polymorphism of benzene and naphthalene with free energy based enhanced molecular dynamics

Prediction and exploration of possible polymorphism in organic crystal compounds are of great importance for industries ranging from organic electronics to pharmaceuticals to high-energy materials. Here we apply our crystal structure prediction procedure and the enhanced molecular dynamics based sampling approach called the Crystal-Adiabatic Free Energy Dynamics (Crystal-AFED) method to benzene and naphthalene. Crystal-AFED allows the free energy landscape of structures to be explored efficiently at any desired temperature and pressure. For each system, we successfully predict the most stable crystal structures at atmospheric pressure and explore the relative Gibbs free energies of predicted polymorphs at high pressures. Using Crystal-AFED sampling, we find that mixed structures, which typically cannot be discovered by standard crystal structure prediction methods, are prevalent in the solid forms of these compounds at high pressure.

Organic crystal polymorphism: a benchmark for dispersion-corrected mean-field electronic structure methods

We analyze the energy landscape of the sixth crystal structure prediction blind test targets with various first principles and semi-empirical quantum chemical methodologies. A new benchmark set of 59 crystal structures (termed POLY59) for testing quantum chemical methods based on the blind test target crystals is presented. We focus on different means to include London dispersion interactions within the density functional theory (DFT) framework. We show the impact of pairwise dispersion corrections like the semi-empirical D2 scheme, the Tkatchenko-Scheffler (TS) method, and the density-dependent dispersion correction dDsC. Recent methodological progress includes higher-order contributions in both the many-body and multipole expansions. We use the D3 correction with Axilrod-Teller-Muto type three-body contribution, the TS based many-body dispersion (MBD), and the nonlocal van der Waals density functional (vdW-DF2). The density functionals with D3 and MBD correction provide an energy ranking of the blind test polymorphs in excellent agreement with the experimentally found structures. As a computationally less demanding method, we test our recently presented minimal basis Hartree-Fock method (HF-3c) and a density functional tight-binding Hamiltonian (DFTB). Considering the speed-up of three to four orders of magnitudes, the energy ranking provided by the low-cost methods is very reasonable. We compare the computed geometries with the corresponding X-ray data where TPSS-D3 performs best. The importance of zero-point vibrational energy and thermal effects on crystal densities is highlighted.

Transferable force field for crystal structure predictions, investigation of performance and exploration of different rescoring strategies using DFT-D methods

A new force field, here called AZ-FF, aimed at being used for crystal structure predictions, has been developed. The force field is transferable to a new type of chemistry without additional training or modifications. This makes the force field very useful in the prediction of crystal structures of new drug molecules since the time-consuming step of developing a new force field for each new molecule is circumvented. The accuracy of the force field was tested on a set of 40 drug-like molecules and found to be very good where observed crystal structures are found at the top of the ranked list of tentative crystal structures. Re-ranking with dispersion-corrected density functional theory (DFT-D) methods further improves the scoring. After DFT-D geometry optimization the observed crystal structure is found at the leading top of the ranking list. DFT-D methods and force field methods have been evaluated for use in predicting properties such as phase transitions upon heating, mechanical properties or intrinsic crystalline solubility. The utility of using crystal structure predictions and the associated material properties in risk assessment in connection with form selection in the drug development process is discussed.

Using crystal structure prediction to rationalize the hydration propensities of substituted adamantane hydrochloride salts

The crystal energy landscapes of the salts of two rigid pharmaceutically active molecules reveal that the experimental structure of amantadine hydrochloride is the most stable structure with the majority of low-energy structures adopting a chain hydrogen-bond motif and packings that do not have solvent accessible voids. By contrast, memantine hydrochloride which differs in the substitution of two methyl groups on the adamantane ring has a crystal energy landscape where all structures within 10 kJ mol−1of the global minimum have solvent-accessible voids ranging from 3 to 14% of the unit-cell volume including the lattice energy minimum that was calculated after removing water from the hydrated memantine hydrochloride salt structure. The success in using crystal structure prediction (CSP) to rationalize the different hydration propensities of these substituted adamantane hydrochloride salts allowed us to extend the model to predict under blind test conditions the experimental crystal structures of the previously uncharacterized 1-(methylamino)adamantane base and its corresponding hydrochloride salt. Although the crystal structure of 1-(methylamino)adamantane was correctly predicted as the second ranked structure on the static lattice energy landscape, the crystallization of aZ′ = 3 structure of 1-(methylamino)adamantane hydrochloride reveals the limits of applying CSP when the contents of the crystallographic asymmetric unit are unknown.