Composite and Low-Cost Approaches for Molecular Crystal Structure Prediction

How Accurate Can Crystal Structure Predictions Be for High-Energy Molecular Crystals?

Molecular crystals have shallow potential energy landscapes, with multiple local minima separated by very small differences in total energy. Predicting molecular packing and molecular conformation in the crystal generally requires ab initio methods of high accuracy, especially when polymorphs are involved. We used dispersion-corrected density functional theory (DFT-D) to assess the capabilities of an evolutionary algorithm (EA) for the crystal structure prediction (CSP) of well-known but challenging high-energy molecular crystals (HMX, RDX, CL-20, and FOX-7). While providing the EA with the experimental conformation of the molecule quickly re-discovers the experimental packing, it is more realistic to start instead from a naïve, flat, or neutral initial conformation, which reflects the limited experimental knowledge we generally have in the computational design of molecular crystals. By doing so, and using fully flexible molecules in fully variable unit cells, we show that the experimental structures can be predicted in fewer than 20 generations. Nonetheless, one must be aware that some molecular crystals have naturally hindered evolutions, requiring as many attempts as there are space groups of interest to predict their structures, and some may require the accuracy of all-electron calculations to discriminate between closely ranked structures. To save resources in this computationally demanding process, we showed that a hybrid xTB/DFT-D approach could be considered in a subsequent study to push the limits of CSP beyond 200+ atoms and for cocrystals.

Discriminating and understanding molecular crystal polymorphism

Polymorph discrimination for a molecular crystal has long been a challenging task, which, nonetheless, is a major concern in the pharmaceutical industry. In this work, we have investigated polymorph discrimination on three different molecular crystals, tetrolic acid, oxalic acid, and oxalyl dihydrazide, covering both packing polymorphism and conformational polymorphism. To gain more understanding, we have performed energy decomposition analysis based on many-body expansion, and have compared the results from the XO-PBC method, that is, the eXtended ONIOM method (XO) with the periodic boundary condition (PBC), with those from some commonly used dispersion corrected density functional theory (DFT-D) methods. It is shown here that, with the XYG3 doubly hybrid functional chosen as the target high level to capture the intra- and short-range intermolecular interactions, and the periodic PBE as the basic low level to take long range interactions into account, the XO-PBC(XYG3:PBE) method not only obtains the correct experimental stability orderings, but also predicts reasonable polymorph energy ranges for all three cases. Our results have demonstrated the usefulness of the present theoretical methods, in particular XO-PBC, while highlighted the importance of a better treatment of different kinds of interactions to be beneficial to polymorph control.

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.

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

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

Revisiting stacking interactions in tetrathiafulvalene and selected derivatives using tight-binding quantum chemical calculations and local coupled-cluster method

The engineering of supramolecular architectures needs accurate descriptions of the intermolecular interactions in crystal structures. Tetrathiafulvalene (TTF) is an effective building block used in the construction of promising functional materials. The parallel packing of the neutral TTF-TTF system was studied previously using the high-level quantum chemical method, advancing it as a valuable model system. The recently developed tight-binding quantum chemical method GFN2-xTB and local coupled-cluster method DLPNO-CCSD(T) were used to investigate the stacking interactions of TTF and selected derivatives deposited in the Cambridge Structural Database. Using the interaction energy of the TTF-TTF dimer calculated at the CCSD(T)/CBS level as the reference, the accuracies of the two methods are investigated. The energy decomposition analysis within the DLPNO-CCSD(T) framework reveals the importance of dispersion interaction in the TTF-related stacking systems. The dispersion interaction density plot vividly shows the magnitude and distribution of the dispersion interaction, providing a revealing insight into the stacking interactions in crystal structures. The results show that the GFN2-xTB and DLPNO-CCSD(T) methods could achieve accuracy at an affordable computational cost, which would be valuable in understanding the nature of parallel stacking in supramolecular systems.

Interplay between London Dispersion, Hubbard U, and Metastable States for Uranium Compounds

High-throughput computational studies of lanthanide and actinide chemistry with density-functional theory are complicated by the need for Hubbard U corrections, which ensure localization of the f-electrons, but can lead to metastable states. This work presents a systematic investigation of the effects of both Hubbard U value and metastable states on the predicted structural and thermodynamic properties of four uranium compounds central to the field of nuclear fuels: UC, UN, UO₂, and UCl₃. We also assess the impact of the exchange-hole dipole moment (XDM) dispersion correction on the computed properties. Overall, the choice of Hubbard U value and inclusion of a dispersion correction cause larger variations in the computed geometric properties than result from metastable states. The weak dependence of structure optimization on metastable states should simplify future high-throughput calculations on actinides. Conversely, addition of the dispersion correction is found to offset the repulsion introduced by the Hubbard U term and provides greatly improved agreement with experiment for both cell volumes and heats of formation. The XDM dispersion correction is largely invariant to the chosen U value, making it a robust dispersion correction for actinide systems.

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

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

Theoretical investigation of amino-acid adsorption on hydroxylated quartz surfaces: dispersion can determine enantioselectivity

Chiral mineral surfaces, such as quartz, are attractive substrates for use in enantioselective separation and may have contributed to the origin of biological homochirality. In this work, we apply density-functional theory and the exchange-hole dipole moment (XDM) dispersion model to study the adsorption of 5 amino acids (glycine, serine, alanine, valine, and phenylalanine) on a hydroxylated α-quartz (0001) surface. It is demonstrated that London dispersion is responsible for 30-50% of the total adsorption energies and its inclusion or omission can reverse predictions of enantioselectivity. Differing dispersion stabilization, caused by the opposing side-chain placements relative to the quartz surface, lead to differences of 1.0 and 1.8 kcal mol^-1 in the adsorption energies of the alanine and phenylalanine enantiomers, respectively. These results are consistent with a 3-point model, with the hydrogen-bonding sites conserved and variations in the dispersion interactions determining enantioselectivity.

First-Principles-Based Force Field for 2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105)

2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105) is a highly promising energetic material (EM) with high safety. Understanding its microscopic response mechanisms within the external stimulus is meaningful for the design of EMs. In order to comprehend the complicated phenomena, it is necessary to employ molecular simulation methods to investigate the response mechanisms with the force field (FF) at an atomic level. In this work, we developed a tailored FF for LLM-105 based on first-principles calculations. The validity of the FF was evaluated by molecular dynamics simulations. The structural parameters of LLM-105 predicted by FF are in good agreement with the experimental values, such as lattice constant, bond length, bond angle, dihedral angle and center of mass, and so forth. Moreover, the FF possesses good performance to describe the structural response on pressure accurately. In general, our work not only builds a balanced FF in gas and condensed phases, but also provides a useful tool to study the properties about LLM-105 at a large scale, which is helpful to improve the understanding about the balance between energy and safety in EMs.

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.

Exploring the applicability of density functional tight binding to transition metal ions. Parameterization for nickel with the spin-polarized DFTB3 model

In this work, we explore the applicability and limitations of the current third order density functional tight binding (DFTB3) formalism for treating transition metal ions using nickel as an example. To be consistent with recent parameterization of DFTB3 for copper, the parametrization for nickel is conducted in a spin-polarized formulation and with orbital-resolved Hubbard parameters and their charge derivatives. The performance of the current parameter set is evaluated based on structural and energetic properties of a set of nickel-containing compounds that involve biologically relevant ligands. Qualitatively similar to findings in previous studies of copper complexes, the DFTB3 results are more reliable for nickel complexes with neutral ligands than for charged ligands; nevertheless, encouraging agreement is noted in comparison to the reference method, B3LYP/aug-cc-pVTZ, especially for structural properties, including cases that exhibit Jahn-Teller distortions; the structures also compare favorably to available X-ray data in the Cambridge Crystallographic Database for a number of nickel-containing compounds. As to limitations, we find it is necessary to use different d shell Hubbard charge derivatives for Ni(I) and Ni(II), due to the distinct electronic configurations for the nickel ion in the respective complexes, and substantial errors are observed for ligand binding energies, especially for charged ligands, d orbital splitting energies and splitting between singlet and triplet spin states for Ni(II) compounds. These observations highlight that future improvement in intra-d correlation and ligand polarization is required to enable the application of the DFTB3 model to complex transition metal ions. © 2018 Wiley Periodicals, Inc.

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