Reliable crystal structure predictions from first principles

Catching the π-Stacks: Prediction of Aggregate Structures of Porphyrin

π-π interactions decisively shape the supramolecular structure and functionality of π-conjugated molecular semiconductor materials. Despite the customizable molecular building blocks, predicting their supramolecular structure remains a challenge. Traditionally, force field methods have been used due to the complexity of these structures, but advances in computational power have enabled ab initio approaches such as density functional theory (DFT). DFT is particularly suitable for finding energetically favorable structures of dye aggregates, which are determined by a large number of different interactions, but a systematic aggregate search can still be very challenging due to the large number of possible geometries. In this work, we show ways to overcome this challenge. We investigate how finely translational and rotational lattices must be structured to identify all energetic minima of π-stack structures, focusing on porphyrins as a prototype challenge. Our approach involves single-point DFT calculations of systematically varied dimer geometries, identification of local energy minima, hierarchical grouping of geometrically similar structures, and optimization of the energetically favorable representatives of each geometric family. This ab initio method provides a general framework for the systematic prediction of aggregate structures and reveals geometrically diverse and energetically favorable dimers.

Equivariant Neural Networks Utilizing Molecular Clusters for Accurate Molecular Crystal Lattice Energy Predictions

Equivariant neural networks have emerged as prominent models in advancing the construction of interatomic potentials due to their remarkable data efficiency and generalization capabilities for out-of-distribution data. Here, we expand the utility of these networks to the prediction of crystal structures consisting of organic molecules. Traditional methods for computing crystal structure properties, such as plane-wave quantum chemical methods based on density functional theory (DFT), are prohibitively resource-intensive, often necessitating compromises in accuracy and the choice of exchange-correlation functional. We present an approach that leverages the efficiency, and transferability of equivariant neural networks, specifically Allegro, to predict molecular crystal structure energies at a reduced computational cost. Our neural network is trained on molecular clusters using a highly accurate Gaussian-type orbital (GTO)-based method as the target level of theory, eliminating the need for costly periodic DFT calculations, while providing access to all families of exchange-corelation functionals and post-Hartree-Fock methods. The trained model exhibits remarkable accuracy in predicting lattice energies, aligning closely with those computed by plane-wave based DFT methods, thus representing significant cost reductions. Furthermore, the Allegro network was seamlessly integrated with the USPEX framework, accelerating the discovery of low-energy crystal structures during crystal structure prediction.

Contributions beyond direct random-phase approximation in the binding energy of solid ethane, ethylene, and acetylene

The random-phase approximation (RPA) includes a subset of higher than second-order correlation-energy contributions, but stays in the same complexity class as the second-order Møller-Plesset perturbation theory (MP2) in both Gaussian-orbital and plane-wave codes. This makes RPA a promising ab initio electronic structure approach for the binding energies of molecular crystals. Still, some issues stand out in practical applications of RPA. Notably, compact clusters of nonpolar molecules are poorly described, and the interaction energies strongly depend on the reference single-determinant state. Using the many-body expansion of the binding energy of a crystal, we investigate those issues and the effect of beyond-RPA corrections. We find the beneficial effect of quartic-scaling exchange and non-ring coupled-cluster doubles corrections. The nonadditive interactions in compact trimers of molecules are improved by using the self-consistent Hartree-Fock orbitals instead of the usual Kohn-Sham states, but this kind of orbital input also leads to underestimated dimer energies. Overall, a substantial improvement over the RPA with a renormalized singles approach is possible at a modest quartic-scaling cost, which encourages further research into additional RPA corrections.

Crystal Polymorph Search in the NPT Ensemble via a Deposition/Sublimation Alchemical Path

The formulation of active pharmaceutical ingredients involves discovering stable crystal packing arrangements or polymorphs, each of which has distinct pharmaceutically relevant properties. Traditional experimental screening techniques utilizing various conditions are commonly supplemented with in silico crystal structure prediction (CSP) to inform the crystallization process and mitigate risk. Predictions are often based on advanced classical force fields or quantum mechanical calculations that model the crystal potential energy landscape but do not fully incorporate temperature, pressure, or solution conditions during the search procedure. This study proposes an innovative alchemical path that utilizes an advanced polarizable atomic multipole force field to predict crystal structures based on direct sampling of the NPT ensemble. The use of alchemical (i.e., nonphysical) intermediates, a novel Monte Carlo barostat, and an orthogonal space tempering bias combine to enhance the sampling efficiency of the deposition/sublimation phase transition. The proposed algorithm was applied to 2-((4-(2-(3,4-dichlorophenyl)ethyl)phenyl)amino)benzoic acid (Cambridge Crystallography Database Centre ID: XAFPAY) as a case study to showcase the algorithm. Each experimentally determined polymorph with one molecule in the asymmetric unit was successfully reproduced via approximately 1000 short 1 ns simulations per space group where each simulation was initiated from random rigid body coordinates and unit cell parameters. Utilizing two threads of a recent Intel CPU (a Xeon Gold 6330 CPU at 2.00 GHz), 1 ns of sampling using the polarizable AMOEBA force field can be acquired in 4 h (equating to more than 300 ns/day using all 112 threads/56 cores of a dual CPU node) within the Force Field X software (https://ffx.biochem.uiowa.edu). These results demonstrate a step forward in the rigorous use of the NPT ensemble during the CSP search process and open the door to future algorithms that incorporate solution conditions using continuum solvation methods.

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.

Ab Initio Crystal Structure Prediction of the Energetic Materials LLM-105, RDX, and HMX

Crystal structure prediction (CSP) is performed for the energetic materials (EMs) LLM-105 and α-RDX, as well as the α and β conformational polymorphs of 1,3,5,7-tetranitro-1,3,5,7-tetraazacyclooctane (HMX), using the genetic algorithm (GA) code, GAtor, and its associated random structure generator, Genarris. Genarris and GAtor successfully generate the experimental structures of all targets. GAtor's symmetric crossover scheme, where the space group symmetries of parent structures are treated as genes inherited by offspring, is found to be particularly effective. However, conducting several GA runs with different settings is still important for achieving diverse samplings of the potential energy surface. For LLM-105 and α-RDX, the experimental structure is ranked as the most stable, with all of the dispersion-inclusive density functional theory (DFT) methods used here. For HMX, the α form was persistently ranked as more stable than the β form, in contrast to experimental observations, even when correcting for vibrational contributions and thermal expansion. This may be attributed to insufficient accuracy of dispersion-inclusive DFT methods or to kinetic effects not considered here. In general, the ranking of some putative structures is found to be sensitive to the choice of the DFT functional and the dispersion method. For LLM-105, GAtor generates a putative structure with a layered packing motif, which is desirable thanks to its correlation with low sensitivity. Our results demonstrate that CSP is a useful tool for studying the ubiquitous polymorphism of EMs and shows promise of becoming an integral part of the EM development pipeline.

Revisiting Homochiral versus Heterochiral Interactions through a Long Detective Story of a Useful Azobis-Nitrile and Puzzling Racemate

This paper documents and reinvestigates the solid-state and crystal structures of 4,4'-azobis-4-cyanopentanoic acid (ACPA), a water-soluble azobis-nitrile of immense utility as a radical initiator in living polymerizations and a labile mechanophore that can be embedded within long polymer chains to undergo selective scission under mechanical activation. Surprisingly, for such applications, both the commercially available reagent and their derivatives are used as "single initiators" when this azonitrile is actually a mixture of stereoisomers. Although the racemate and meso compounds were identified more than half a century ago and their enantiomers were separated by classical resolution, there have been confusing narratives dealing with their characterization, the existence of a conglomeratic phase, and fractional crystallization. Our results report on the X-ray crystal structures of all stereoisomers for the first time, along with further details on enantiodiscrimination and the always intriguing arguments accounting for the stability of homochiral versus heterochiral crystal aggregates. To this end, metadynamic (MTD) simulations on stereoisomer molecular aggregates were performed to capture the incipient nucleation events at the picosecond time scale. This analysis sheds light on the driving homochiral aggregation of ACPA enantiomers.

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.

Efficient Crystal Structure Prediction for Structurally Related Molecules with Accurate and Transferable Tailor-Made Force Fields

Crystal structure prediction (CSP) his generally used to complement experimental solid form screening and applied to individual molecules in drug development. The fast development of algorithms and computing resources offers the opportunity to use CSP earlier and for a broader range of applications in the drug design cycle. This study presents a novel paradigm of CSP specifically designed for structurally related molecules, referred to as Quick-CSP. The approach prioritizes more accurate physics through robust and transferable tailor-made force fields (TMFFs), such that significant efficiency gains are achieved through the reduction of expensive ab initio calculations. The accuracy of the TMFF is increased by the introduction of electrostatic multipoles, and the fragment-based force field parameterization scheme is demonstrated to be transferable for a family of chemically related molecules. The protocol is benchmarked with structurally related compounds from the Bromodomain and Extraterminal (BET) domain inhibitors series. A new convergence criterion is introduced that aims at performing only as many ab initio optimizations of crystal structures as required to locate the bottom of the crystal energy landscape within a user-defined accuracy. The overall approach provides significant cost savings ranging from three- to eight-fold less than the full-CSP workflow. The reported advancements expand the scope and utility of the underlying CSP building blocks as well as their novel reassembly to other applications earlier in the drug design cycle to guide molecule design and selection.