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

Reducing the computational cost of NMR crystallography of organic powders at natural isotopic abundance with the help of 13 C-13 C dipolar couplings

Structure determination of functional organic compounds remains a formidable challenge when the sample exists as a powder. Nuclear magnetic resonance crystallography approaches based on the comparison of experimental and Density Functional Theory (DFT)-computed ¹ H chemical shifts have already demonstrated great potential for structure determination of organic powders, but limitations still persist. In this study, we discuss the possibility of using ¹³ C-¹³ C dipolar couplings quantified on powdered theophylline at natural isotopic abundance with the help of dynamic nuclear polarization, to realize a DFT-free, rapid screening of a pool of structures predicted by ab initio random structure search. We show that although ¹³ C-¹³ C dipolar couplings can identify structures possessing long range structural motifs and unit cell parameters close to those of the true structure, it must be complemented with other data to recover information about the presence and the chemical nature of the supramolecular interactions.

Completing the picture of tolazamide polymorphism under extreme conditions: a low-temperature study

We present the results of an experimental and computational study of structural changes in two polymorphs of tolazamide {systematic name: 1-[(azepan-1-ylamino)carbonyl]-4-methylbenzenesulfonamide}, C₁₄H₂₁N₃O₃S, on cooling to 100 K and reverse heating. No phase transitions occurred in this temperature range. The anisotropy of the thermal expansion was different for the two polymorphs and differed from that reported previously for the hydrostatic compression. The changes in different intermolecular contacts responsible for the strain anisotropy were analysed. Relative shortening of the contacts was related directly to their initial length and reversely to the steric density around them. Increasing steric density is likely to be the driving force for the conformational ordering of the azepane ring under compression.

Nature of halogen-centered intermolecular interactions in crystal growth and design: Fluorine-centered interactions in dimers in crystalline hexafluoropropylene as a prototype

The wide occurrence of halogen-centered noncovalent interactions in crystal growth and design prompted this study, which includes a mini review of recent advances in the field. Particular emphasis is placed on providing compelling theoretical evidence of the formation of these interactions between sites of positive electrostatic potential, as well as between sites of negative electrostatic potential, localized on the electrostatic surfaces of the bound fluorine atoms in a prototypical system, hexafluoropropylene (C₃ F₆ ), upon its interaction with another same molecule to form (C₃ F₆ )₂ dimers. The existence of σ- and π-hole interactions is shown for the stable dimers. Even so, weakly bound interactions locally responsible in holding the molecular fragments together cannot and should not be overlooked since they are partly responsible for determining the overall geometry of the crystal. The results of combined quantum theory of atoms in molecules, molecular electrostatic surface potential, and reduced density gradient noncovalent interaction analyses showed that these latter interactions do indeed play a role in the stability and growth of crystalline C₃ F₆ itself and the (C₃ F₆ )₂ dimers. A symmetry adapted perturbation theory energy decomposition analysis leads to the conclusion that a great majority of the (C₃ F₆ )₂ dimers examined are the consequence of dispersion (and electrostatics), with nonnegligible contribution from polarization, which together competes with an exchange repulsion component to determine the equilibrium geometries. In a few structures of the (C₃ F₆ )₂ dimer, the fluorine is found to serve as a six-center five-bond donor/acceptor, as found for carbon in other systems (Malischewski and Seppelt, Angew. Chem. Int. Ed. 2017, 56, 368). © 2019 Wiley Periodicals, Inc.

A full additive QM/MM scheme for the computation of molecular crystals with extension to many-body expansions

An additive quantum mechanics/molecular mechanics (QM/MM) model for the theoretical investigation of molecular crystals (AC-QM/MM) is presented. At the one-body level, a single molecule is chosen as the QM region. The MM region around it consists of a finite cluster of explicit MM atoms, represented by point charges and Lennard-Jones potentials, with additional background charges to mimic periodic electrostatics. Cluster charges are QM-derived and calculated self-consistently to ensure a polarizable embedding. We have also considered the extension to many-body QM corrections, calculating the interactions of a central molecule to neighboring units in the crystal. Full gradient expressions have been derived, also including symmetry information. The scheme allows for the calculation of molecular properties as well as unconstrained optimizations of the molecular geometry and cell parameters with respect to the lattice energy. Benchmarking the approach with the X23 reference set confirms the convergence pattern of the many-body extension although a comparison to plane-wave density functional theory reveals a systematic overestimation of cohesive energies by 6-16 kJ mol^-1. While the scheme primarily aims to provide an inexpensive and flexible way to model a molecule in a crystal environment, it can also be used to reach highly accurate cohesive energies by the straightforward application of wave function correlated approaches. Calculations with local coupled cluster with singles, doubles, and perturbative triples, albeit limited to numerical gradients, show an impressive agreement with experimental estimates for small molecular crystals.

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.

Soft Crystals: Flexible Response Systems with High Structural Order

A new material concept of soft crystals is proposed. Soft crystals respond to gentle stimuli such as vapor exposure and rubbing but maintain their structural order and exhibit remarkable visual changes in their shape, color, and luminescence. Various interesting examples of soft crystals are introduced in the article. By exploring the interesting formation and phase-transition phenomena of soft crystals through interdisciplinary collaboration, new materials having both the characteristics of ordered hard crystals and those of flexible soft matter are expected.

Machine-Learned Fragment-Based Energies for Crystal Structure Prediction

Crystal structure prediction involves a search of a complex configurational space for local minima corresponding to stable crystal structures, which can be performed efficiently using atom-atom force fields for the assessment of intermolecular interactions. However, for challenging systems, the limitations in the accuracy of force fields prevent a reliable assessment of the relative thermodynamic stability of potential structures, while the cost of fully quantum mechanical approaches can limit applications of the methods. We present a method to rapidly improve force field lattice energies by correcting two-body interactions with a higher level of theory in a fragment-based approach and predicting these corrections with machine learning. Corrected lattice energies with commonly used density functionals and second order perturbation theory (MP2) all significantly improve the ranking of experimentally known polymorphs where the rigid molecule model is applicable. The relative lattice energies of known polymorphs are also found to systematically improve with the fragment corrections. Predicting two-body interactions with atom-centered symmetry functions in a Gaussian process is found to give highly accurate results using as little as 10-20% of the data for training, reducing the cost of the energy correction by up to an order of magnitude. The machine learning approach opens up the possibility of more widespread use of fragment-based methods in crystal structure prediction, whose increased accuracy at a low computational cost will benefit applications in areas such as polymorph screening and computer-guided materials discovery.

Challenges in the use of atomistic simulations to predict solubilities of drug-like molecules

Background: Solubility is a physical property of high importance to the pharmaceutical industry, the prediction of which for potential drugs has so far been a hard task. We attempted to predict the solubility of acetylsalicylic acid (ASA) by estimating the absolute chemical potentials of its most stable polymorph and of solutions with different concentrations of the drug molecule. Methods: Chemical potentials were estimated from all-atom molecular dynamics simulations.  We used the Einstein molecule method (EMM) to predict the absolute chemical potential of the solid and solvation free energy calculations to predict the excess chemical potentials of the liquid-phase systems. Results: Reliable estimations of the chemical potentials for the solid and for a single ASA molecule using the EMM required an extremely large number of intermediate states for the free energy calculations, meaning that the calculations were extremely demanding computationally. Despite the computational cost, however, the computed value did not agree well with the experimental value, potentially due to limitations with the underlying energy model. Perhaps better values could be obtained with a better energy model; however, it seems likely computational cost may remain a limiting factor for use of this particular approach to solubility estimation.    Conclusions: Solubility prediction of drug-like solids remains computationally challenging, and it appears that both the underlying energy model and the computational approach applied may need improvement before the approach is suitable for routine use.

Understanding non-covalent interactions in larger molecular complexes from first principles

Non-covalent interactions pervade all matter and play a fundamental role in layered materials, biological systems, and large molecular complexes. Despite this, our accumulated understanding of non-covalent interactions to date has been mainly developed in the tens-of-atoms molecular regime. This falls considerably short of the scales at which we would like to understand energy trends, structural properties, and temperature dependencies in materials where non-covalent interactions have an appreciable role. However, as more reference information is obtained beyond moderately sized molecular systems, our understanding is improving and we stand to gain pertinent insights by tackling more complex systems, such as supramolecular complexes, molecular crystals, and other soft materials. In addition, accurate reference information is needed to provide the drive for extending the predictive power of more efficient workhorse methods, such as density functional approximations that also approximate van der Waals dispersion interactions. In this perspective, we discuss the first-principles approaches that have been used to obtain reference interaction energies for beyond modestly sized molecular complexes. The methods include quantum Monte Carlo, symmetry-adapted perturbation theory, non-canonical coupled cluster theory, and approaches based on the random-phase approximation. By considering the approximations that underpin each method, the most accurate theoretical references for supramolecular complexes and molecular crystals to date are ascertained. With these, we also assess a handful of widely used exchange-correlation functionals in density functional theory. The discussion culminates in a framework for putting into perspective the accuracy of high-level wavefunction-based methods and identifying future challenges.

Polymorphism in p-aminobenzoic acid

We review the polymorphism of p-aminobenzoic acid (pABA), a model drug compound whose crystallisation and polymorphic behaviour has been extensively studied in recent years.

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.

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.

Explaining crystallization preferences of two polyphenolic diastereoisomers by crystal structure prediction

Crystal structure prediction is used to understand the differences in crystallization of catechin and epicatechin, and to explore the predictability of solvate formation.

The interplay among molecular structures, crystal symmetries and lattice energy landscapes revealed using unsupervised machine learning: a closer look at pyrrole azaphenacenes

Using unsupervised machine learning and CSPs to help crystallographers better understand how crystallizations are affected by molecular structures.

Reliable and practical computational description of molecular crystal polymorphs

Reliable prediction of the polymorphic energy landscape of a molecular crystal would yield profound insight into drug development in terms of the existence and likelihood of late-appearing polymorphs. However, the computational prediction of molecular crystal polymorphs is highly challenging due to the high dimensionality of conformational and crystallographic space accompanied by the need for relative free energies to within 1 kJ/mol per molecule. In this study, we combine the most successful crystal structure sampling strategy with the most successful first-principles energy ranking strategy of the latest blind test of organic crystal structure prediction methods. Specifically, we present a hierarchical energy ranking approach intended for the refinement of relative stabilities in the final stage of a crystal structure prediction procedure. Such a combined approach provides excellent stability rankings for all studied systems and can be applied to molecular crystals of pharmaceutical importance.

The role of molecular modelling and simulation in the discovery and deployment of metal-organic frameworks for gas storage and separation

Metal-organic frameworks (MOFs) are highly tuneable, extended-network, crystalline, nanoporous materials with applications in gas storage, separations, and sensing. We review how molecular models and simulations of gas adsorption in MOFs have informed the discovery of performant MOFs for methane, hydrogen, and oxygen storage, xenon, carbon dioxide, and chemical warfare agent capture, and xylene enrichment. Particularly, we highlight how large, open databases of MOF crystal structures, post-processed to enable molecular simulations, are a platform for computational materials discovery. We discuss how to orient research efforts to routinise the computational discovery of MOFs for adsorption-based engineering applications.

A jumping crystal predicted with molecular dynamics and analysed with TLS refinement against powder diffraction data

By running a temperature series of molecular dynamics (MD) simulations starting from the known low-temperature phase, the experimentally observed phase transition in a 'jumping crystal' was captured, thereby providing a prediction of the unknown crystal structure of the high-temperature phase and clarifying the phase-transition mechanism. The phase transition is accompanied by a discontinuity in two of the unit-cell parameters. The structure of the high-temperature phase is very similar to that of the low-temperature phase. The anisotropic displacement parameters calculated from the MD simulations readily identified libration as the driving force behind the phase transition. Both the predicted crystal structure and the phase-transition mechanism were verified experimentally using TLS (translation, libration, screw) refinement against X-ray powder diffraction data.

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.

Metastable silica high pressure polymorphs as structural proxies of deep Earth silicate melts

Modelling of processes involving deep Earth liquids requires information on their structures and compression mechanisms. However, knowledge of the local structures of silicates and silica (SiO2) melts at deep mantle conditions and of their densification mechanisms is still limited. Here we report the synthesis and characterization of metastable high-pressure silica phases, coesite-IV and coesite-V, using in situ single-crystal X-ray diffraction and ab initio simulations. Their crystal structures are drastically different from any previously considered models, but explain well features of pair-distribution functions of highly densified silica glass and molten basalt at high pressure. Built of four, five-, and six-coordinated silicon, coesite-IV and coesite-V contain SiO6 octahedra, which, at odds with 3rd Pauling's rule, are connected through common faces. Our results suggest that possible silicate liquids in Earth's lower mantle may have complex structures making them more compressible than previously supposed.

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.

Repulsion-dispersion parameters for the modelling of organic molecular crystals containing N, O, S and Cl

In lattice energy models that combine ab initio and empirical components, it is important to ensure consistency between these components so that meaningful quantitative results are obtained. A method for deriving parameters of atom-atom repulsion dispersion potentials for crystals, tailored to different ab initio models, is presented. It is based on minimization of the sum of squared deviations between experimental and calculated structures and energies. The solution algorithm is designed to avoid convergence to local minima in the parameter space by combining a deterministic low-discrepancy sequence for the generation of multiple initial parameter guesses with an efficient local minimization algorithm. The proposed approach is applied to derive transferable exp-6 potential parameters suitable for use in conjunction with a distributed multipole electrostatics model derived from isolated molecule charge densities calculated at the M06/6-31G(d,p) level of theory. Data for hydrocarbons, azahydrocarbons, oxohydrocarbons, organosulphur compounds and chlorohydrocarbons are used for the estimation. A good fit is achieved for the new set of parameters with a mean absolute error in sublimation enthalpies of 4.1 kJ mol-1 and an average rmsd15 of 0.31 Å. The parameters are found to perform well on a separate cross-validation set of 39 compounds.

Pervasive Delocalisation Error Causes Spurious Proton Transfer in Organic Acid-Base Co-Crystals

Dispersion-corrected density-functional theory (DFT-D) methods have become the workhorse of many computational protocols for molecular crystal structure prediction due to their efficiency and convenience. However, certain limitations of DFT, such as delocalisation error, are often overlooked or are too expensive to remedy in solid-state applications. This error can lead to artificial stabilisation of charge transfer and, in this work, it is found to affect the correct identification of the protonation site in multicomponent acid-base crystals. As such, commonly used DFT-D methods cannot be applied with any reliability to the study of acid-base co-crystals or salts, while hybrid functionals remain too restrictive for routine use. This presents an impetus for the development of new functionals with reduced delocalisation error for solid-state applications; the structures studied herein constitute an excellent benchmark for this purpose.

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.

NMR Crystallography of the Polymorphs of Metergoline

Two polymorphs of the drug compound metergoline (C25H29N3O2) were investigated in detail by solid-state NMR measurements. The results have been analysed by an advanced procedure, which uses experimental input together with the results of quantum chemical calculations that were performed for molecular crystals. In this way, it was possible to assign the total of 40 1H–13C correlation pairs in a highly complex system, namely, in the dynamically disordered polymorph with two independent molecules in the unit cell of a large volume of 4234 Å3. For the simpler polymorph, which exhibits only small-amplitude motions and has just one molecule in the unit cell with a volume of 529.0 Å3, the values of the principal elements of the 13C chemical shift tensors were measured. Additionally, for this polymorph, a set of crystal structure predictions were generated, and the {13C, 1H} isotropic and 13C anisotropic chemical shielding data were computed while using the gauge-including projector augmented-wave approach combined with the “revised Perdew-Burke-Ernzerhof“ exchange-correlation functional (GIPAW-RPBE). The experimental and theoretical results were combined in an application of the newly developed strategy to polymorph discrimination. This research thus opens up new routes towards more accurate characterization of the polymorphism of drug formulations.

Preconditioners for the geometry optimisation and saddle point search of molecular systems

A class of preconditioners is introduced to enhance geometry optimisation and transition state search of molecular systems. We start from the Hessian of molecular mechanical terms, decompose it and retain only its positive definite part to construct a sparse preconditioner matrix. The construction requires only the computation of the gradient of the corresponding molecular mechanical terms that are already available in popular force field software packages. For molecular crystals, the preconditioner can be combined straightforwardly with the exponential preconditioner recently introduced for periodic systems. The efficiency is demonstrated on several systems using empirical, semiempirical and ab initio potential energy surfaces.

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.

Azulene revisited: solid-state structure, invariom modeling and lattice-energy minimization of a classical example of disorder

The molecular and solid-state structure of azulene both raise fundamental questions. Therefore, the disordered crystal structure of azulene was re-refined with invariom non-spherical atomic scattering factors from new single-crystal X-ray diffraction data with a resolution of d = 0.45 Å. An unconstrained refinement results in a molecular geometry with C s symmetry. Refinements constrained to fulfill C 2v symmetry, as observed in the gas phase and in high-level ab initio calculations, lead to similar figures of merit and residual densities as unconstrained ones. Such models are consistent with the structures from microwave spectroscopy and electron diffraction, albeit they are not the same. It is shown that for the disorder present in azulene, the invariom model describes valence electron density as successfully as it does for non-disordered structures, although the disorder still leads to high correlations mainly between positional parameters. Lattice-energy minimizations on a variety of ordered model structures using dispersion-corrected DFT calculations reveal that the local deviations from the average structure are small. Despite the molecular dipole moment there is no significant molecular ordering in any spatial direction. A superposition of all ordered model structures leads to a calculated average structure, which explains not only the experimental determined atomic coordinates, but also the apparently unusual experimental anisotropic displacement parameters.

Tailoring van der Waals dispersion interactions with external electric charges

van der Waals (vdW) dispersion interactions strongly impact the properties of molecules and materials. Often, the description of vdW interactions should account for the coupling with pervasive electric fields, stemming from membranes, ionic channels, liquids, or nearby charged functional groups. However, this quantum-mechanical effect has been omitted in atomistic simulations, even in widely employed electronic-structure methods. Here, we develop a model and study the effects of an external charge on long-range vdW correlations. We show that a positive external charge stabilizes dispersion interactions, whereas a negative charge has an opposite effect. Our analytical results are benchmarked on a series of (bio)molecular dimers and supported by calculations with high-level correlated quantum-chemical methods, which estimate the induced dispersion to reach up to 35% of intermolecular binding energy (4 kT for amino-acid dimers at room temperature). Our analysis bridges electrostatic and electrodynamic descriptions of intermolecular interactions and may have implications for non-covalent reactions, exfoliation, dissolution, and permeation through biological membranes.

Genarris: Random generation of molecular crystal structures and fast screening with a Harris approximation

We present Genarris, a Python package that performs configuration space screening for molecular crystals of rigid molecules by random sampling with physical constraints. For fast energy evaluations, Genarris employs a Harris approximation, whereby the total density of a molecular crystal is constructed via superposition of single molecule densities. Dispersion-inclusive density functional theory is then used for the Harris density without performing a self-consistency cycle. Genarris uses machine learning for clustering, based on a relative coordinate descriptor developed specifically for molecular crystals, which is shown to be robust in identifying packing motif similarity. In addition to random structure generation, Genarris offers three workflows based on different sequences of successive clustering and selection steps: the "Rigorous" workflow is an exhaustive exploration of the potential energy landscape, the "Energy" workflow produces a set of low energy structures, and the "Diverse" workflow produces a maximally diverse set of structures. The latter is recommended for generating initial populations for genetic algorithms. Here, the implementation of Genarris is reported and its application is demonstrated for three test cases.

DNP-enhanced solid-state NMR spectroscopy of active pharmaceutical ingredients

Solid-state NMR spectroscopy has become a valuable tool for the characterization of both pure and formulated active pharmaceutical ingredients (APIs). However, NMR generally suffers from poor sensitivity that often restricts NMR experiments to nuclei with favorable properties, concentrated samples, and acquisition of one-dimensional (1D) NMR spectra. Here, we review how dynamic nuclear polarization (DNP) can be applied to routinely enhance the sensitivity of solid-state NMR experiments by one to two orders of magnitude for both pure and formulated APIs. Sample preparation protocols for relayed DNP experiments and experiments on directly doped APIs are detailed. Numerical spin diffusion models illustrate the dependence of relayed DNP enhancements on the relaxation properties and particle size of the solids and can be used for particle size determination when the other factors are known. We then describe the advanced solid-state NMR experiments that have been enabled by DNP and how they provide unique insight into the molecular and macroscopic structure of APIs. For example, with large sensitivity gains provided by DNP, natural isotopic abundance, ¹³ C-¹³ C double-quantum single-quantum homonuclear correlation NMR spectra of pure APIs can be routinely acquired. DNP also enables solid-state NMR experiments with unreceptive quadrupolar nuclei such as ² H, ¹⁴ N, and ³⁵ Cl that are commonly found in APIs. Applications of DNP-enhanced solid-state NMR spectroscopy for the molecular level characterization of low API load formulations such as commercial tablets and amorphous solid dispersions are described. Future perspectives for DNP-enhanced solid-state NMR experiments on APIs are briefly discussed.

Molecular dynamics based enhanced sampling of collective variables with very large time steps

Enhanced sampling techniques that target a set of collective variables and that use molecular dynamics as the driving engine have seen widespread application in the computational molecular sciences as a means to explore the free-energy landscapes of complex systems. The use of molecular dynamics as the fundamental driver of the sampling requires the introduction of a time step whose magnitude is limited by the fastest motions in a system. While standard multiple time-stepping methods allow larger time steps to be employed for the slower and computationally more expensive forces, the maximum achievable increase in time step is limited by resonance phenomena, which inextricably couple fast and slow motions. Recently, we introduced deterministic and stochastic resonance-free multiple time step algorithms for molecular dynamics that solve this resonance problem and allow ten- to twenty-fold gains in the large time step compared to standard multiple time step algorithms [P. Minary et al., Phys. Rev. Lett. 93, 150201 (2004); B. Leimkuhler et al., Mol. Phys. 111, 3579-3594 (2013)]. These methods are based on the imposition of isokinetic constraints that couple the physical system to Nosé-Hoover chains or Nosé-Hoover Langevin schemes. In this paper, we show how to adapt these methods for collective variable-based enhanced sampling techniques, specifically adiabatic free-energy dynamics/temperature-accelerated molecular dynamics, unified free-energy dynamics, and by extension, metadynamics, thus allowing simulations employing these methods to employ similarly very large time steps. The combination of resonance-free multiple time step integrators with free-energy-based enhanced sampling significantly improves the efficiency of conformational exploration.

First principles crystal engineering of nonlinear optical materials. I. Prototypical case of urea

The crystalline materials with nonlinear optical (NLO) properties are critically important for several technological applications, including nanophotonic and second harmonic generation devices. Urea is often considered to be a standard NLO material, due to the combination of non-centrosymmetric crystal packing and capacity for intramolecular charge transfer. Various approaches to crystal engineering of non-centrosymmetric molecular materials were reported in the literature. Here we propose using global lattice energy minimization to predict the crystal packing from the first principles. We developed a methodology that includes the following: (1) parameter derivation for polarizable force field AMOEBA; (2) local minimizations of crystal structures with these parameters, combined with the evolutionary algorithm for a global minimum search, implemented in program USPEX; (3) filtering out duplicate polymorphs produced; (4) reoptimization and final ranking based on density functional theory (DFT) with many-body dispersion (MBD) correction; and (5) prediction of the second-order susceptibility tensor by finite field approach. This methodology was applied to predict virtual urea polymorphs. After filtering based on packing similarity, only two distinct packing modes were predicted: one experimental and one hypothetical. DFT + MBD ranking established non-centrosymmetric crystal packing as the global minimum, in agreement with the experiment. Finite field approach was used to predict nonlinear susceptibility, and H-bonding was found to account for a 2.5-fold increase in molecular hyperpolarizability to the bulk value.

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.

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.

Grain boundary phases in bcc metals

We report a computational discovery of novel grain boundary structures and multiple grain boundary phases in elemental body-centered cubic (bcc) metals represented by tungsten, tantalum and molybdenum. While grain boundary structures created by the γ-surface method as a union of two perfect half crystals have been studied extensively, it is known that the method has limitations and does not always predict the correct ground states. Herein, we use a newly developed computational tool, based on evolutionary algorithms, to perform a grand-canonical search of high-angle symmetric tilt and twist boundaries, and we find new ground states and multiple phases that cannot be described using the conventional structural unit model. We use molecular dynamics (MD) simulations to demonstrate that the new structures can coexist at finite temperature in a closed system, confirming that these are examples of different grain boundary phases. The new ground state is confirmed by first-principles calculations.

Evaluation of General and Tailor Made Force Fields via X-ray Thermal Diffuse Scattering Using Molecular Dynamics and Monte Carlo Simulations of Crystalline Aspirin

We have performed a comparison of the experimental thermal diffuse scattering (TDS) from crystalline Aspirin (form I) to that calculated from molecular dynamics (MD) simulations based on a variety of general force fields and a tailor-made force field (TMFF). A comparison is also made with Monte Carlo (MC) simulations which use a "harmonic network" approach to describe the intermolecular interactions. These comparisons were based on the hypothesis that TDS could be a useful experimental data in validation of such simulation parameter sets, especially when calculations of dynamical properties (e.g., thermodynamic free energies) from molecular crystals are concerned. Currently such a validation of force field parameters against experimental data is often limited to calculation of specific physical properties, e.g., absolute lattice energies usually at 0 K or heat capacity measurements. TDS harvested from in-house or synchrotron experiments comprises highly detailed structural information representative of the dynamical motions of the crystal lattice. Thus, TDS is a well-suited experimental data-driven means of cross validating theoretical approaches targeted at understanding dynamical properties of crystals. We found from the results of our investigation that the TMFF and COMPASS (from the commercial software "Materials Studio") parameter sets gave the best agreement with experiment. From our homologous MC simulation analysis we are able to show that force constants associated with the molecular torsion angles are likely to be a strong contributing factor for the apparent reason why these aforementioned force fields performed better.

Computational prediction of drug solubility in water-based systems: Qualitative and quantitative approaches used in the current drug discovery and development setting

In this review we will discuss recent advances in computational prediction of solubility in water-based solvents. Our focus is set on recent advances in predictions of biorelevant solubility in media mimicking the human intestinal fluids and on new methods to predict the thermodynamic cycle rather than prediction of solubility in pure water through quantitative structure property relationships (QSPR). While the literature is rich in QSPR models for both solubility and melting point, a physicochemical property strongly linked to the solubility, recent advances in the modelling of these properties make use of theory and computational simulations to better predict these properties or processes involved therein (e.g. solid state crystal lattice packing, dissociation of molecules from the lattice and solvation). This review serves to provide an update on these new approaches and how they can be used to more accurately predict solubility, and also importantly, inform us on molecular interactions and processes occurring during drug dissolution and solubilisation.

Composite and Low-Cost Approaches for Molecular Crystal Structure Prediction

Molecular crystal structure prediction (CSP) requires evaluating differences in lattice energy between candidate crystal structures accurately and efficiently. In this work, we explore and compare several low-cost alternatives to dispersion-corrected density-functional theory (DFT) in the plane-waves/pseudopotential approximation, the most accurate and general approach used for CSP at present. Three types of low-cost methods are considered: DFT with a small basis set of finite-support numerical orbitals (the SIESTA method), dispersion-corrected Gaussian small or minimal-basis-set Hartree-Fock and DFT with additional empirical corrections (HF-3c and PBEh-3c), and self-consistent-charge dispersion-corrected density-functional tight binding (SCC-DFTB3-D3). In addition, we study the performance of composite methods that comprise a geometry optimization using a low-cost approach followed by a single-point calculation using the accurate but comparatively expensive B86bPBE-XDM functional. All methods were tested for their abilities to produce absolute lattice energies, relative lattice energies, and crystal geometries. We show that assessing various methods by their ability to produce absolute lattice energies can be misleading and that relative lattice energies are a much better indicator of performance in CSP. The EE14 set of relative solubilities of homochiral and heterochiral chiral crystals is proposed for relative lattice-energy benchmarking. Our results show that PBE-D2 plus a DZP basis set of numerical orbitals coupled with a final B86bPBE-XDM single-point calculation gives excellent performance at a fraction of the cost of a full B86bPBE-XDM calculation, although the results are sensitive to the particular details of the computational protocol. The B86bPBE-XDM//PBE-D2/DZP method was subsequently tested in a practical CSP application from our recent work on the crystal structure of the enantiopure and racemate forms of 1-aza[6]helicene, a chiral organic semiconductor. Our results show that this multilevel method is able to correctly reproduce the energy ranking of both crystal forms.

GAtor: A First-Principles Genetic Algorithm for Molecular Crystal Structure Prediction

We present the implementation of GAtor, a massively parallel, first-principles genetic algorithm (GA) for molecular crystal structure prediction. GAtor is written in Python and currently interfaces with the FHI-aims code to perform local optimizations and energy evaluations using dispersion-inclusive density functional theory (DFT). GAtor offers a variety of fitness evaluation, selection, crossover, and mutation schemes. Breeding operators designed specifically for molecular crystals provide a balance between exploration and exploitation. Evolutionary niching is implemented in GAtor by using machine learning to cluster the dynamically updated population by structural similarity and then employing a cluster-based fitness function. Evolutionary niching promotes uniform sampling of the potential energy surface by evolving several subpopulations, which helps overcome initial pool biases and selection biases (genetic drift). The various settings offered by GAtor increase the likelihood of locating numerous low-energy minima, including those located in disconnected, hard to reach regions of the potential energy landscape. The best structures generated are re-relaxed and re-ranked using a hierarchy of increasingly accurate DFT functionals and dispersion methods. GAtor is applied to a chemically diverse set of four past blind test targets, characterized by different types of intermolecular interactions. The experimentally observed structures and other low-energy structures are found for all four targets. In particular, for Target II, 5-cyano-3-hydroxythiophene, the top ranked putative crystal structure is a Z' = 2 structure with P1̅ symmetry and a scaffold packing motif, which has not been reported previously.