Accurate Modeling of Organic Molecular Crystals by Dispersion-Corrected Density Functional Tight Binding (DFTB)

Computational predictions of cocrystal formation: A benchmark study of 28 assemblies comparing five methods from high-throughput to advanced models

Cocrystals are assemblies of more than one type of molecule stabilized through noncovalent interactions. They are promising materials for improved drug formulation in which the stability, solubility, or biocompatibility of the active pharmaceutical ingredient (API) is improved by including a coformer. In this work, a range of density functional theory (DFT) and density functional tight binding (DFTB) models are systematically compared for their ability to predict the lattice enthalpy of a broad range of existing pharmaceutically relevant cocrystals. These range from cocrystals containing model compounds 4,4'-bipyridine and oxalic acid to those with the well benchmarked APIs of aspirin and paracetamol, all tested with a large set of alternative coformers. For simple cocrystals, there is a general consensus in lattice enthalpy calculated by the different DFT models. For the cocrystals with API coformers the cocrystals, enthalpy predictions depend strongly on the DFT model. The significantly lighter DFTB models predict unrealistic values of lattice enthalpy even for simple cocrystals.

A Molecular View of Methane Activation on Ni(111) through Enhanced Sampling and Machine Learning

A combination of machine learned interatomic potentials (MLIPs) and enhanced sampling simulations is used to investigate the activation of methane on a Ni(111) surface. The work entails the development and iterative refinement of MLIPs, initially trained on a data set constructed via ab initio molecular dynamics simulations, supplemented by adaptive biasing forces, to enrich the sampling of catalytically relevant configurations. Our results reveal that upon incorporation of collective variables that capture the behavior of the reactant molecule, as well as additional frames that describe the dynamic response of the catalytic surface, it is possible to enhance considerably the accuracy of predicted energies and forces. By employing enhanced sampling schemes in the refinement of the MLIP, we systematically explore the potential energy surface, leading to a refined MLIP capable of predicting density functional theory-level energies and forces and replicating key geometric characteristics of the catalytic system. The resulting free energy landscapes at several temperatures provide a detailed view of the thermodynamics and dynamics of methane activation. Specifically, as methane approaches and dissociates on the catalytic surface, the process involves the dynamic interplay of CH4 and the Ni catalyst that includes both enthalpic and entropic contributions. The progression toward the transition state involves a CH4 moiety that is increasingly restrained in its ability to rotate or translate, while the stage following the transition state is characterized by a notable rise of the Ni atom that interacts with the cleaved C-H bond. This leads to an increase in the mobility of the adsorbed species, a feature that becomes more pronounced at higher temperatures.

Condensed Matter Systems Exposed to Radiation: Multiscale Theory, Simulations, and Experiment

This roadmap reviews the new, highly interdisciplinary research field studying the behavior of condensed matter systems exposed to radiation. The Review highlights several recent advances in the field and provides a roadmap for the development of the field over the next decade. Condensed matter systems exposed to radiation can be inorganic, organic, or biological, finite or infinite, composed of different molecular species or materials, exist in different phases, and operate under different thermodynamic conditions. Many of the key phenomena related to the behavior of irradiated systems are very similar and can be understood based on the same fundamental theoretical principles and computational approaches. The multiscale nature of such phenomena requires the quantitative description of the radiation-induced effects occurring at different spatial and temporal scales, ranging from the atomic to the macroscopic, and the interlinks between such descriptions. The multiscale nature of the effects and the similarity of their manifestation in systems of different origins necessarily bring together different disciplines, such as physics, chemistry, biology, materials science, nanoscience, and biomedical research, demonstrating the numerous interlinks and commonalities between them. This research field is highly relevant to many novel and emerging technologies and medical applications.

Computer-Aided Design and Biological Evaluation of Diazaspirocyclic D4R Antagonists

Parkinson's disease (PD) is a neurodegenerative disorder characterized by the progressive loss of dopaminergic neurons in the substantia nigra, resulting in motor dysfunction. Current treatments are primarily centered around enhancing dopamine signaling or providing dopamine replacement therapy and face limitations such as reduced efficacy over time and adverse side effects. To address these challenges, we identified selective dopamine receptor subtype 4 (D4R) antagonists not previously reported as potential adjuvants for PD management. In this study, a library screening and artificial neural network quantitative structure-activity relationship (QSAR) modeling with experimentally driven library design resulted in a class of spirocyclic compounds to identify candidate D4R antagonists. However, developing selective D4R antagonists suitable for clinical translation remains a challenge.

PyRETIS 3: Conquering rare and slow events without boundaries

We present and discuss the advancements made in PyRETIS 3, the third instalment of our Python library for an efficient and user-friendly rare event simulation, focused to execute molecular simulations with replica exchange transition interface sampling (RETIS) and its variations. Apart from a general rewiring of the internal code towards a more modular structure, several recently developed sampling strategies have been implemented. These include recently developed Monte Carlo moves to increase path decorrelation and convergence rate, and new ensemble definitions to handle the challenges of long-lived metastable states and transitions with unbounded reactant and product states. Additionally, the post-analysis software PyVisa is now embedded in the main code, allowing fast use of machine-learning algorithms for clustering and visualising collective variables in the simulation data.

Exploring the Mechanisms behind Non-aromatic Fluorescence with the Density Functional Tight Binding Method

Recent experimental findings reveal nonconventional fluorescence emission in biological systems devoid of conjugated bonds or aromatic compounds, termed non-aromatic fluorescence (NAF). This phenomenon is exclusive to aggregated or solid states and remains absent in monomeric solutions. Previous studies focused on small model systems in vacuum show that the carbonyl stretching mode along with strong interaction of short hydrogen bonds (SHBs) remains the primary vibrational mode explaining NAF in these systems. In order to simulate larger model systems taking into account the effects of the surrounding environment, in this work we propose using the density functional tight-binding (DFTB) method in combination with non-adiabatic molecular dynamics (NAMD) and the mixed quantum/molecular mechanics (QM/MM) approach. We investigate the mechanism behind NAF in the crystal structure of l-pyroglutamine-ammonium, comparing it with the related nonfluorescent amino acid l-glutamine. Our results extend our previous findings to more realistic systems, demonstrating the efficiency and robustness of the proposed DFTB method in the context of NAMD in biological systems. Furthermore, due to its inherent low computational cost, this method allows for a better sampling of the nonradiative events at the conical intersection which is crucial for a complete understanding of this phenomenon. Beyond contributing to the ongoing exploration of NAF, this work paves the way for future application of this method in more complex biological systems such as amyloid aggregates, biomaterials, and non-aromatic proteins.

Effect of interlayer slipping on the geometric, thermal and adsorption properties of 2D covalent organic frameworks: a comprehensive review based on computational modelling studies

Two-dimensional covalent organic frameworks (2D-COFs) are a class of crystalline porous organic polymers, consisting of 2D-planar sheets stacked together perpendicularly via noncovalent forces. Since their discovery, 2D-COFs have attracted extensive attention for optoelectronic and adsorption applications. Owing to the layer stacking nature of 2D COFs, various new slipped structures that are energetically favourable can be designed. These interlayer slipped structures are actively responsible for tuning (mostly enhancing) the optoelectronic properties, thermal properties, and mechanical strength of 2D COFs. This review summarizes the effect of interlayer slipping on the energetic stability, electronic behaviour and gas adsorption properties of 2D layered COFs, which is explained through computational modelling simulations. Since computational modelling offers a deep insight into electronic behaviour at the atomic scale, which is potentially impossible through experimental techniques, the introduction and role of computational techniques in such studies have also been described.

From Organic Fragments to Photoswitchable Catalysts: The OFF-ON Structural Repository for Transferable Kernel-Based Potentials

Structurally and conformationally diverse databases are needed to train accurate neural networks or kernel-based potentials capable of exploring the complex free energy landscape of flexible functional organic molecules. Curating such databases for species beyond "simple" drug-like compounds or molecules composed of well-defined building blocks (e.g., peptides) is challenging as it requires thorough chemical space mapping and evaluation of both chemical and conformational diversities. Here, we introduce the OFF-ON (organic fragments from organocatalysts that are non-modular) database, a repository of 7869 equilibrium and 67,457 nonequilibrium geometries of organic compounds and dimers aimed at describing conformationally flexible functional organic molecules, with an emphasis on photoswitchable organocatalysts. The relevance of this database is then demonstrated by training a local kernel regression model on a low-cost semiempirical baseline and comparing it with a PBE0-D3 reference for several known catalysts, notably the free energy surfaces of exemplary photoswitchable organocatalysts. Our results demonstrate that the OFF-ON data set offers reliable predictions for simulating the conformational behavior of virtually any (photoswitchable) organocatalyst or organic compound composed of H, C, N, O, F, and S atoms, thereby opening a computationally feasible route to explore complex free energy surfaces in order to rationalize and predict catalytic behavior.

Mechanoresponsive Flexible Crystals

Flexible crystals have gained significant attention owing to their remarkable pliability, plasticity, and adaptability, making them highly popular in various research and application fields. The main challenges in developing flexible crystals lie in the rational design, preparation, and performance optimization of such crystals. Therefore, a comprehensive understanding of the fundamental origins of crystal flexibility is crucial for establishing evaluation criteria and design principles. This Perspective offers a retrospective analysis of the development of flexible crystals over the past two decades. It summarizes the elastic standards and possible plastic bending mechanisms tailored to diverse flexible crystals and analyzes the assessment of their theoretical basis and applicability. Meanwhile, the compatibility between crystal elasticity and plasticity has been discussed, unveiling the immense prospects of elastic/plastic crystals for applications in biomedicine, flexible electronic devices, and flexible optics. Furthermore, this Perspective presents state-of-the-art experimental avenues and analysis methods for investigating molecular interactions in molecular crystals, which is vital for the future exploration of the mechanisms of crystal flexibility.

Addressing electronic and dynamical evolution of molecules and molecular clusters: DFTB simulations of energy relaxation in polycyclic aromatic hydrocarbons

We present a review of the capabilities of the density functional based Tight Binding (DFTB) scheme to address the electronic relaxation and dynamical evolution of molecules and molecular clusters following energy deposition via either collision or photoabsorption. The basics and extensions of DFTB for addressing these systems and in particular their electronic states and their dynamical evolution are reviewed. Applications to PAH molecules and clusters, carbonaceous systems of major interest in astrochemical/astrophysical context, are reported. A variety of processes are examined and discussed such as collisional hydrogenation, fast collisional processes and induced electronic and charge dynamics, collision-induced fragmentation, photo-induced fragmentation, relaxation in high electronic states, electronic-to-vibrational energy conversion and statistical versus non-statistical fragmentation. This review illustrates how simulations may help to unravel different relaxation mechanisms depending on various factors such as the system size, specific electronic structure or excitation conditions, in close connection with experiments.

Accelerated Organic Crystal Structure Prediction with Genetic Algorithms and Machine Learning

We present a high-throughput, end-to-end pipeline for organic crystal structure prediction (CSP)─the problem of identifying the stable crystal structures that will form from a given molecule based only on its molecular composition. Our tool uses neural network potentials to allow for efficient screening and structural relaxation of generated crystal candidates. Our pipeline consists of two distinct stages: random search, whereby crystal candidates are randomly generated and screened, and optimization, where a genetic algorithm (GA) optimizes this screened population. We assess the performance of each stage of our pipeline on 21 molecules taken from the Cambridge Crystallographic Data Centre's CSP blind tests. We show that random search alone yields matches for ≈50% of targets. We then validate the potential of our full pipeline, making use of the GA to optimize the root-mean-square deviation between crystal candidates and the experimentally derived structure. With this approach, we are able to find matches for ≈80% of candidates with 10-100 times smaller initial population sizes than when using random search. Lastly, we run our full pipeline with an ANI model that is trained on a small data set of molecules extracted from crystal structures in the Cambridge Structural Database, generating ≈60% of targets. By leveraging machine learning models trained to predict energies at the density functional theory level, our pipeline has the potential to approach the accuracy of ab initio methods and the efficiency of empirical force fields.

Benchmark Data Set of Crystalline Organic Semiconductors

This work reports a Benchmark Data set of Crystalline Organic Semiconductors to test calculations of the structural and electronic properties of these materials in the solid state. The data set contains 67 crystals consisting of mostly rigid molecules with a single dominant conformer, covering the majority of known structural types. The experimental crystal structure is available for the entire data set, whereas zero-temperature unit cell volume can be reliably estimated for a subset of 28 crystals. Using this subset, we benchmark r2SCAN-D3 and PBE-D3 density functionals. Then, for the entire data set, we benchmark approximate density functional theory (DFT) methods, including GFN1-xTB and DFTB3(3ob-3-1), with various dispersion corrections against r2SCAN-D3. Our results show that r2SCAN-D3 geometries are accurate within a few percent, which is comparable to the statistical uncertainty of experimental data at a fixed temperature, but the unit cell volume is systematically underestimated by 2% on average. The several times faster PBE-D3 provides an unbiased estimate of the volume for all systems except for molecules with highly polar bonds, for which the volume is substantially overestimated in correlation with the underestimation of atomic charges. Considered approximate DFT methods are orders of magnitude faster and provide qualitatively correct but overcompressed crystal structures unless the dispersion corrections are fitted by unit cell volume.

Development of a Crystallization-Induced Diastereomer Transformation of Oxime Isomers for the Asymmetric Synthesis of (1S,6R)-3,9-Diazabicyclo[4.2.1]nonane

Herein we report a practical crystallization-induced diastereomer transformation (CIDT) of oxime isomers for the scalable asymmetric synthesis of the bicyclic diamine (1S,6R)-3,9-diazabicyclo[4.2.1]nonane derivative that serves as a valuable building block in medicinal chemistry. The developed approach utilizes (S)-phenylethylamine as a chiral auxiliary handle for CIDT, and the starting nortropinone derivative is prepared in one step from commercially available materials. The resulting E-oxime is subjected to a stereospecific Beckmann rearrangement, followed by reduction of the resulting lactam with LiAlH4 to afford the monoprotected (1S,6R)-3,9-diazabicyclo[4.2.1]nonane derivative. The development of the CIDT and understanding of the mechanistic implications leading to the high selectivity are reported.

Two excited-state datasets for quantum chemical UV-vis spectra of organic molecules

We present two open-source datasets that provide time-dependent density-functional tight-binding (TD-DFTB) electronic excitation spectra of organic molecules. These datasets represent predictions of UV-vis absorption spectra performed on optimized geometries of the molecules in their electronic ground state. The GDB-9-Ex dataset contains a subset of 96,766 organic molecules from the original open-source GDB-9 dataset. The ORNL_AISD-Ex dataset consists of 10,502,904 organic molecules that contain between 5 and 71 non-hydrogen atoms. The data reveals the close correlation between the magnitude of the gaps between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), and the excitation energy of the lowest singlet excited state energies quantitatively. The chemical variability of the large number of molecules was examined with a topological fingerprint estimation based on extended-connectivity fingerprints (ECFPs) followed by uniform manifold approximation and projection (UMAP) for dimension reduction. Both datasets were generated using the DFTB+ software on the "Andes" cluster of the Oak Ridge Leadership Computing Facility (OLCF).

Weak noncovalent interactions in two positional isomers of acrylonitrile derivatives: inputs from PIXEL energy, Hirshfeld surface and QTAIM analyses

A single crystal X-ray diffraction analysis was performed on two positional isomers (m-tolyl and p-tolyl) of acrylonitrile derivatives, namely, (Z)-3-(4-(pyridin-2-yl) phenyl)-2-(m-tolyl) acrylonitrile (1) and (Z)-3-(4-(pyridin-2-yl)phenyl)-2-(p-tolyl) acrylonitrile (2). Compound 1 crystallized in the monoclinic P2₁/n space group with two crystallographically independent molecules. Compound 2 also possesses two crystallographically independent molecules and crystallized in the triclinic P-1 space group. The Hirshfeld surface analysis revealed that, in both isomers, intermolecular H⋅⋅⋅H/C/N contacts contribute significantly to the crystal packing. More than 40% of the contribution arises from intermolecular C-H⋅⋅⋅C(π) contacts. In both compounds, the relative contribution of these contacts is comparable, indicating that the positional isomeric effects are marginal. The structures in which these isomers are arranged in the solid state are very similar, and the lattice energies are also comparable between the isomers. The Coulomb-London-Pauli-PIXEL (CLP-PIXEL) energy analysis identified the energetically significant dimers. The strength of the intra- and intermolecular interactions was evaluated using the quantum theory of atoms in molecules approach. The UV-Vis absorbance in three different solvents (chloroform, ethanol, and ethyl acetate) for isomers 1 and 2 are very similar. This result is in good agreement with the time-dependent density-functional theory (TD-DFT) calculations.

Adaptive Peptide Molecule as the Promising Highly-Efficient Gas-Sensor Material: In Silico Study

Gas sensors are currently employed in various applications in fields such as medicine, ecology, and food processing, and serve as monitoring tools for the protection of human health, safety, and quality of life. Herein, we discuss a promising direction in the research and development of gas sensors based on peptides-biomolecules with high selectivity and sensitivity to various gases. Thanks to the technique developed in this work, which uses a framework based on the density-functional tight-binding theory (DFTB), the most probable adsorption centers were identified and used to describe the interaction of some analyte molecules with peptides. The DFTB method revealed that the physical adsorption of acetone, ammonium, benzene, ethanol, hexane, methanol, toluene, and trinitrotoluene had a binding energy in the range from -0.28 eV to -1.46 eV. It was found that peptides may adapt to the approaching analyte by changing their volume up to a maximum value of approx. 13%, in order to confine electron clouds around the adsorbed molecule. Based on the results obtained, the prospects for using the proposed peptide configurations in gas sensor devices are good.

Reparameterization of the chemical-potential equalization model with DFTB3: A practical balance between accuracy and transferability

To improve the performance of the third-order density-functional tight-binding method (DFTB3) for non-covalent interactions involving organic and biological molecules, a chemical-potential equalization (CPE) approach was introduced [J. Phys. Chem. A, 116, 9131 (2012)] and parameterized for the H, C, N, O, and S chemical elements [J. Chem. Phys., 143, 084123 (2015)]. Based largely on equilibrium structures, the parameterized DFTB3/CPE models were shown to exhibit improvements in molecular polarizabilities and intermolecular interactions. With more extensive analyses, however, we observe here that the available DFTB3/CPE models have two critical limitations: (1) they lead to sharply varying potential energy surfaces, thus causing numerical instability in molecular dynamics (MD) simulations, and (2) they lead to spurious interactions at short distances for some dimer complexes. These shortcomings are attributed to the employed screening functions and the overfitting of CPE parameters. In this work, we introduce a new strategy to simplify the parameterization procedure and significantly reduce free parameters down to four global (i.e., independent of element type) ones. With this strategy, two new models, DFTB3/CPE(r) and DFTB3/CPE(r†) are parameterized. The new models lead to smooth potential energy surfaces, stable MD simulations, and alleviate the spurious interactions at short distances, thus representing consistent improvements for both neutral and ionic hydrogen bonds.

Surrogate Based Genetic Algorithm Method for Efficient Identification of Low-Energy Peptide Structures

Identification of the most stable structure(s) of a system is a prerequisite for the calculation of any of its properties from first-principles. However, even for relatively small molecules, exhaustive explorations of the potential energy surface (PES) are severely hampered by the dimensionality bottleneck. In this work, we address the challenging task of efficiently sampling realistic low-lying peptide coordinates by resorting to a surrogate based genetic algorithm (GA)/density functional theory (DFT) approach (sGADFT) in which promising candidates provided by the GA are ultimately optimized with DFT. We provide a benchmark of several computational methods (GAFF, AMOEBApro13, PM6, PM7, DFTB3-D3(BJ)) as possible prescanning surrogates and apply sGADFT to two test case systems that are (i) two isomer families of the protonated Gly-Pro-Gly-Gly tetrapeptide (Masson, A.; J. Am. Soc. Mass Spectrom.2015, 26, 1444-1454) and (ii) the doubly protonated cyclic decapeptide gramicidin S (Nagornova, N. S.; J. Am. Chem. Soc.2010, 132, 4040-4041). We show that our GA procedure can correctly identify low-energy minima in as little as a few hours. Subsequent refinement of surrogate low-energy structures within a given energy threshold (≤10 kcal/mol (i), ≤5 kcal/mol (ii)) via DFT relaxation invariably led to the identification of the most stable structures as determined from high-resolution infrared (IR) spectroscopy at low temperature. The sGADFT method therefore constitutes a highly efficient route for the screening of realistic low-lying peptide structures in the gas phase as needed for instance for the interpretation and assignment of experimental IR spectra.

Regularized by Physics: Graph Neural Network Parametrized Potentials for the Description of Intermolecular Interactions

Simulations of molecular systems using electronic structure methods are still not feasible for many systems of biological importance. As a result, empirical methods such as force fields (FF) have become an established tool for the simulation of large and complex molecular systems. The parametrization of FF is, however, time-consuming and has traditionally been based on experimental data. Recent years have therefore seen increasing efforts to automatize FF parametrization or to replace FF with machine-learning (ML) based potentials. Here, we propose an alternative strategy to parametrize FF, which makes use of ML and gradient-descent based optimization while retaining a functional form founded in physics. Using a predefined functional form is shown to enable interpretability, robustness, and efficient simulations of large systems over long time scales. To demonstrate the strength of the proposed method, a fixed-charge and a polarizable model are trained on ab initio potential-energy surfaces. Given only information about the constituting elements, the molecular topology, and reference potential energies, the models successfully learn to assign atom types and corresponding FF parameters from scratch. The resulting models and parameters are validated on a wide range of experimentally and computationally derived properties of systems including dimers, pure liquids, and molecular crystals.

Importance of Nuclear Quantum Effects on Aqueous Electrolyte Transport under Confinement in Ti3C2 MXenes

Protons display a high chemical activity and strongly affect the charge storage capability in confined interlayer spaces of two-dimensional (2D) materials. As such, an accurate representation of proton dynamics under confinement is important for understanding and predicting charge storage dynamics in these materials. While often ignored in atomistic-scale simulations, nuclear quantum effects (NQEs), e.g., tunneling, can be significant under confinement even at room temperature. Using the thermostatted ring polymer molecular dynamics implementation of path integral molecular dynamics (PIMD) in conjunction with the ReaxFF force field, density functional tight binding (DFTB), and NequIP neural network potential simulations, we investigate the role of NQEs on proton and water transport in bulk water and aqueous electrolytes under confinement in Ti₃C₂ MXenes. Although overall NQEs are relatively small, especially in bulk, we find that they can alter both quantitative values and qualitative trends on both proton transport and water self-diffusion under confinement relative to classical MD predictions. Therefore, our results suggest the need for NQEs to be considered to simulate aqueous systems under confinement for both qualitative and quantitative accuracy.

Molecular Mechanisms of Phosphoester Bond Formation in Water Using Tight-Binding Ab Initio Molecular Dynamics

Phosphate groups are ubiquitous in biomolecules and are usually incorporated through phosphoester bonds between alcohol groups and orthophosphate. The formation of this bond is exceptionally difficult, with associated barriers of 30-45 kcal/mol in the absence of catalysts. In abiotic conditions, polymerizing nucleic acids without enzymes remains very challenging and is still a partly unsolved problem that severely questions the RNA World hypothesis for the origins of life. Offering a solution to this problem would involve a detailed knowledge of the reaction energetics and mechanisms, yet these remain not fully understood at a molecular level, especially because of the very slow reaction rates that represent a significant challenge for the experiments. The number of involved reaction coordinates and the possible role of the solvent in assisting the reaction are challenging for computational studies. Here, we use extensive ab initio molecular dynamics simulations using semiempirical tight-binding methods and enhanced sampling to address these issues. We first show that the choice of the tight-binding method is greatly limited by the instability of the water liquid phase for most DFTB generations and parameter sets that are widely available. We then focus on a model reaction involving methanol and orthophosphate, for which the two protonation states (mono- and dianionic) that are dominant around neutral pH are considered. We compare different proton coordinates that enable (or not) the participation of solvent water molecules. Our simulations suggest that in all cases, a dissociative associative mechanism, with an intermediate metaphosphate, is favored. The main difference between the two phosphate species is that reaction with the monoanion is assisted by the substrate, while that with the dianion involves solvent water molecules. Our results are in agreement with early experimental measurements, but the reaction barriers are underestimated in our framework. We believe that our approach provides an interesting perspective on how to sample the reaction phase space efficiently, but it calls for future studies using more accurate descriptions of chemical reactivity.

Thermal decomposition mechanisms of energetic CL-20-based co-crystals: quantum molecular dynamics simulations

The decomposition mechanisms of energetic CL-20:2,4-dinitro-2,4-diazapentane (DNP) and CL-20:2,4-dinitro-2,4-diazaheptane (DNG) co-crystals at high temperatures (1000, 2000, and 3000 K) were studied by density functional tight-binding molecular dynamics (DFTB-MD) simulations. At different temperatures, their decomposition mechanisms are very different. At 1000 K, conformational changes are observed only for the CL-20:DNG co-crystal, in which the CL-20 changes from β-CL-20 to γ-CL-20. When the temperature is increased to 2000 K, CL-20, DNP, and DNG begin to decompose, and there are five paths for the main initial mechanisms. Further increasing the temperature to 3000 K promotes a more complete decomposition. The initial reactions of CL-20 in the two co-crystals have two channels. There are two initial decomposition channels in the DNP molecule and only one channel in the DNG molecule. As the temperature increases, the decomposition products of the two co-crystals are different. Our work may provide the in-depth understanding of the decomposition mechanisms of high-energy CL-20-based co-crystals at high temperatures.

Computational methods for calculation of protein-ligand binding affinities in structure-based drug design

Drug design is an expensive and time-consuming process. Any method that allows reducing the time the costs of the drug development project can have great practical value for the pharmaceutical industry. In structure-based drug design, affinity prediction methods are of great importance. The majority of methods used to predict binding free energy in protein-ligand complexes use molecular mechanics methods. However, many limitations of these methods in describing interactions exist. An attempt to go beyond these limits is the application of quantum-mechanical description for all or only part of the analyzed system. However, the extensive use of quantum mechanical (QM) approaches in drug discovery is still a demanding challenge. This chapter briefly reviews selected methods used to calculate protein-ligand binding affinity applied in virtual screening (VS), rescoring of docked poses, and lead optimization stage, including QM methods based on molecular simulations.

Quantum chemical investigation of the effect of alkali metal ions on the dynamic structure of water in aqueous solutions

We report quantum chemical molecular dynamics (MD) simulations based on the density-functional tight-binding (DFTB) method to investigate the effect of K+, Na+, and Mg2+ ions in aqueous solutions on the static and dynamic structure of bulk water at room temperature and with various concentrations. The DFTB/MD simulations were validated for the description of ion solvation in aqueous ionic solutions by comparing static pair distribution functions (PDFs) as well as the cation solvation shell between experimental and available ab initio DFT data. The effect of the cations on the water structure, as well as relative differences between K+, Na+, and Mg2+ cations, were analyzed in terms of atomically resolved PDFs as well as time-dependent Van Hove correlation functions (VHFs). The investigation of the VHFs reveals that salt ions generally slow down the dynamic decay of the pair correlations in the water solvation sphere, irrespective of the cation size or charge. The analysis of partial metal-oxygen VHFs indicates that there are long-lived correlations between water and Na+ over long distances, in contrast to K+ and Mg2+.

Applications of Artificial Intelligence and Machine Learning Algorithms to Crystallization

Artificial intelligence and specifically machine learning applications are nowadays used in a variety of scientific applications and cutting-edge technologies, where they have a transformative impact. Such an assembly of statistical and linear algebra methods making use of large data sets is becoming more and more integrated into chemistry and crystallization research workflows. This review aims to present, for the first time, a holistic overview of machine learning and cheminformatics applications as a novel, powerful means to accelerate the discovery of new crystal structures, predict key properties of organic crystalline materials, simulate, understand, and control the dynamics of complex crystallization process systems, as well as contribute to high throughput automation of chemical process development involving crystalline materials. We critically review the advances in these new, rapidly emerging research areas, raising awareness in issues such as the bridging of machine learning models with first-principles mechanistic models, data set size, structure, and quality, as well as the selection of appropriate descriptors. At the same time, we propose future research at the interface of applied mathematics, chemistry, and crystallography. Overall, this review aims to increase the adoption of such methods and tools by chemists and scientists across industry and academia.

Controlling the Formation of Two Concomitant Polymorphs in Hg(II) Coordination Polymers

Controlling the formation of the desired product in the appropriate crystalline form is the fundamental breakthrough of crystal engineering. On that basis, the preferential formation between polymorphic forms, which are referred to as different assemblies achieved by changing the disposition or arrangement of the forming units within the crystalline structure, is one of the most challenging topics still to be understood. Herein, we have observed the formation of two concomitant polymorphs with general formula {[Hg(Pip)₂(4,4′-bipy)]·DMF}_n (P1A, P1B; Pip = piperonylic acid; 4,4′-bipy = 4,4′-bipyridine). Besides, [Hg(Pip)₂(4,4′-bipy)]_n (2) has been achieved during the attempts to isolate these polymorphs. The selective synthesis of P1A and P1B has been successfully achieved by changing the synthetic conditions. The formation of each polymorphic form has been ensured by unit cell measurements and decomposition temperature. The elucidation of their crystal structure revealed P1A and P1B as polymorphs, which originates from the Hg(II) cores and intermolecular associations, especially pinpointed by Hg···π and π···π interactions. Density functional theory (DFT) calculations suggest that P1B, which shows Hg(II) geometries that are further from ideality, is more stable than P1A by 13 kJ·mol^–1 per [Hg(Pip)₂(4,4′-bipy)]·DMF formula unit, and this larger stability of P1B arises mainly from metal···π and π···π interactions between chains. As a result, these structural modifications lead to significant variations of their solid-state photoluminescence.

We have successfully isolated two concomitant polymorphs with formula {[Hg(Pip)₂(4,4′-bipy)]·DMF}_n (P1A and P1B), as well as their desolvated form 2. Then, both polymorphs were selectively synthesized by temperature or anion-template formation. Their crystal structures revealed distorted pentagonal pyramidal geometries and show that differences arise from geometry and packing that led to different solid-state photoluminescence emissions. According to periodic-DFT calculations, distortions in P1B are counterbalanced leading to a more stable form by Hg(II)···π and π···π interactions.

A complete description of thermodynamic stabilities of molecular crystals

Predicting stable polymorphs of molecular crystals remains one of the grand challenges of computational science. Current methods invoke approximations to electronic structure and statistical mechanics and thus fail to consistently reproduce the delicate balance of physical effects determining thermodynamic stability. We compute the rigorous ab initio Gibbs free energies for competing polymorphs of paradigmatic compounds, using machine learning to mitigate costs. The accurate description of electronic structure and full treatment of quantum statistical mechanics allow us to predict the experimentally observed phase behavior. This constitutes a key step toward the first-principles design of functional materials for applications from photovoltaics to pharmaceuticals.

Predictions of relative stabilities of (competing) molecular crystals are of great technological relevance, most notably for the pharmaceutical industry. However, they present a long-standing challenge for modeling, as often minuscule free energy differences are sensitively affected by the description of electronic structure, the statistical mechanics of the nuclei and the cell, and thermal expansion. The importance of these effects has been individually established, but rigorous free energy calculations for general molecular compounds, which simultaneously account for all effects, have hitherto not been computationally viable. Here we present an efficient “end to end” framework that seamlessly combines state-of-the art electronic structure calculations, machine-learning potentials, and advanced free energy methods to calculate ab initio Gibbs free energies for general organic molecular materials. The facile generation of machine-learning potentials for a diverse set of polymorphic compounds—benzene, glycine, and succinic acid—and predictions of thermodynamic stabilities in qualitative and quantitative agreement with experiments highlight that predictive thermodynamic studies of industrially relevant molecular materials are no longer a daunting task.

BENCHMARKS OF THE DENSITY FUNCTIONAL TIGHT-BINDING METHOD FOR REDOX, PROTONATION AND ELECTRONIC PROPERTIES OF QUINONES

Organic materials with controllable molecular design and sustainable resources are promising electrode materials. Crystalline quinones have been investigated in a variety of rechargeable battery chemistries due to their ubiquitous nature,...

Structural stability and initial decomposition mechanisms of BTF crystal induced by vacancy defects: a computational study

Density functional tight binding (DFTB) and DFTB-based molecular dynamics (DFTB-MD) were used to study the effects of vacancy defects on the structure, stability, and initial decomposition mechanisms of condensed phase...

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

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

Understanding the aqueous chemistry of quinoline and the diazanaphthalenes: insight from DFT study

The inter-fragment interactions at various binding sites and the overall cluster stability of quinolone (QNOL), cinnoline (CNOL), quinazoline (QNAZ), and quinoxaline (QNOX) complexes with H2O were studied using the density functional theory (DFT) approach. The adsorption and H-bond binding energies, and the energy decomposition mechanism was considered to determine the relative stabilization status of the studied clusters. Scanning tunneling microscopy (STM), natural bonding orbitals (NBO) and charge decomposition were studied to expose the electronic distribution and interaction between fragments. The feasibility of formations of the various complexes were also studied by considering their thermodynamic properties. Results from adsorption studies confirmed the actual adsorption of H2O molecules on the various binding sites studied, with QNOX clusters exhibiting the best adsorptions. Charge decomposition analysis (CDA) revealed significant charge transfer from substrate to H2O fragment in most complexes, except in QNOL, CNOL and QNAZ clusters with H2O at binding position 4, where much charges are back-donated to substrate. The O---H inter-fragment bonds was discovered to be stronger than counterpart N---H bonds in the complexes, whilst polarity indices confirmed N---H as more polar covalent than O---H bonds. Thermodynamic considerations revealed that the formation process of all studied complexes are endothermic (+ve ΔH f ) and non-spontaneous (+ve ΔG f ).

First principles reactive simulation for equation of state prediction

The high cost of density functional theory (DFT) has hitherto limited the ab initio prediction of the equation of state (EOS). In this article, we employ a combination of large scale computing, advanced simulation techniques, and smart data science strategies to provide an unprecedented ab initio performance analysis of the high explosive pentaerythritol tetranitrate (PETN). Comparison to both experiment and thermochemical predictions reveals important quantitative limitations of DFT for EOS prediction and thus the assessment of high explosives. In particular, we find that DFT predicts the energy of PETN detonation products to be systematically too high relative to the unreacted neat crystalline material, resulting in an underprediction of the detonation velocity, pressure, and temperature at the Chapman-Jouguet state. The energetic bias can be partially accounted for by high-level electronic structure calculations of the product molecules. We also demonstrate a modeling strategy for mapping chemical composition across a wide parameter space with limited numerical data, the results of which suggest additional molecular species to consider in thermochemical modeling.

Cis-Trans Isomerization and Thermal Decomposition Mechanisms of a Series of Nx (x = 4, 8, 10, 11) Chain-Catenated Energetic Crystals

Nitrogen-rich compounds based on heteroaromatic rings with different lengths of nitrogen chains are at the forefront of the energetic materials field. We studied the decomposition processes and reaction kinetics of a series of N_x (x = 4, 8, 10, 11) chain-catenated energetic crystals at various temperatures (2400-3000 K) based on a combinational strategy based on density functional tight binding molecular dynamics (DFTB-MD) simulations and density functional theory (DFT). The results show that the thermal decomposition and reaction kinetics are dependent on both the temperature and nitrogen chain's length. There are two sequential stages in the initial decomposition process for the crystals N₈ and N₁₀: (i) competition between cis-trans isomerization and initial unimolecular decomposition and (ii) subsequent complicated global decomposition reactions. Increasing either the temperature or nitrogen chain's length will accelerate the competition and make initial decomposition dominate. However, cis-trans isomerization does not occur in the crystals N₄ and N₁₁. The dominant initiation paths for N₄, N₈, and N₁₀ occur in the heterocycle and in the bond between the heterocycle and azo group, while that for N₁₁ is ring elimination. The decomposition reactions exhibit a clear first-order kinetics character. The energy paths based on DFT calculations are determined as an addition to the DFTB-MD results. Our findings provide insights into the comprehensive understanding of thermal decomposition behaviors of nitrogen chain-catenated and even all-nitrogen energetic materials.

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

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

Machine learning-accelerated quantum mechanics-based atomistic simulations for industrial applications

Atomistic simulations have become an invaluable tool for industrial applications ranging from the optimization of protein-ligand interactions for drug discovery to the design of new materials for energy applications. Here we review recent advances in the use of machine learning (ML) methods for accelerated simulations based on a quantum mechanical (QM) description of the system. We show how recent progress in ML methods has dramatically extended the applicability range of conventional QM-based simulations, allowing to calculate industrially relevant properties with enhanced accuracy, at reduced computational cost, and for length and time scales that would have otherwise not been accessible. We illustrate the benefits of ML-accelerated atomistic simulations for industrial R&D processes by showcasing relevant applications from two very different areas, drug discovery (pharmaceuticals) and energy materials. Writing from the perspective of both a molecular and a materials modeling scientist, this review aims to provide a unified picture of the impact of ML-accelerated atomistic simulations on the pharmaceutical, chemical, and materials industries and gives an outlook on the exciting opportunities that could emerge in the future.

Crystal Structure Prediction Methods for Organic Molecules: State of the Art

The prediction of the crystal structures that a given organic molecule is likely to form is an important theoretical problem of significant interest for the pharmaceutical and agrochemical industries, among others. As evidenced by a series of six blind tests organized over the past 2 decades, methodologies for crystal structure prediction (CSP) have witnessed substantial progress and have now reached a stage of development where they can begin to be applied to systems of practical significance. This article reviews the state of the art in general-purpose methodologies for CSP, placing them within a common framework that highlights both their similarities and their differences. The review discusses specific areas that constitute the main focus of current research efforts toward improving the reliability and widening applicability of these methodologies, and offers some perspectives for the evolution of this technology over the next decade.

A quantitative evaluation of computational methods to accelerate the study of alloxazine-derived electroactive compounds for energy storage

Alloxazines are a promising class of organic electroactive compounds for application in aqueous redox flow batteries (ARFBs), whose redox properties need to be tuned further for higher performance. High-throughput computational screening (HTCS) enables rational and time-efficient study of energy storage compounds. We compared the performance of computational chemistry methods, including the force field based molecular mechanics, semi-empirical quantum mechanics, density functional tight binding, and density functional theory, on the basis of their accuracy and computational cost in predicting the redox potentials of alloxazines. Various energy-based descriptors, including the redox reaction energies and the frontier orbital energies of the reactant and product molecules, were considered. We found that the lowest unoccupied molecular orbital (LUMO) energy of the reactant molecules is the best performing chemical descriptor for alloxazines, which is in contrast to other classes of energy storage compounds, such as quinones that we reported earlier. Notably, we present a flexible in silico approach to accelerate both the singly and the HTCS studies, therewithal considering the level of accuracy versus measured electrochemical data, which is readily applicable for the discovery of alloxazine-derived organic compounds for energy storage in ARFBs.

Data-efficient machine learning for molecular crystal structure prediction

The combination of modern machine learning (ML) approaches with high-quality data from quantum mechanical (QM) calculations can yield models with an unrivalled accuracy/cost ratio. However, such methods are ultimately limited by the computational effort required to produce the reference data. In particular, reference calculations for periodic systems with many atoms can become prohibitively expensive for higher levels of theory. This trade-off is critical in the context of organic crystal structure prediction (CSP). Here, a data-efficient ML approach would be highly desirable, since screening a huge space of possible polymorphs in a narrow energy range requires the assessment of a large number of trial structures with high accuracy. In this contribution, we present tailored Δ-ML models that allow screening a wide range of crystal candidates while adequately describing the subtle interplay between intermolecular interactions such as H-bonding and many-body dispersion effects. This is achieved by enhancing a physics-based description of long-range interactions at the density functional tight binding (DFTB) level-for which an efficient implementation is available-with a short-range ML model trained on high-quality first-principles reference data. The presented workflow is broadly applicable to different molecular materials, without the need for a single periodic calculation at the reference level of theory. We show that this even allows the use of wavefunction methods in CSP.

Structure, intermolecular interactions, and dynamic properties of NTO crystals with impurity defects: a computational study

Frontier orbitals distribute in the position of impurity molecules, whose adjacent NTO molecules begin to decompose first.

Density functional theory predictions of the mechanical properties of crystalline materials

The DFT-predicted mechanical properties of crystalline materials are crucial knowledge for their screening, design, and exploitation.

Comparison of computational chemistry methods for the discovery of quinone-based electroactive compounds for energy storage

High-throughput computational screening (HTCS) is a powerful approach for the rational and time-efficient design of electroactive compounds. The effectiveness of HTCS is dependent on accuracy and speed at which the performance descriptors can be estimated for possibly millions of candidate compounds. Here, a systematic evaluation of computational methods, including force field (FF), semi-empirical quantum mechanics (SEQM), density functional based tight binding (DFTB), and density functional theory (DFT), is performed on the basis of their accuracy in predicting the redox potentials of redox-active organic compounds. Geometry optimizations at low-level theories followed by single point energy (SPE) DFT calculations that include an implicit solvation model are found to offer equipollent accuracy as the high-level DFT methods, albeit at significantly lower computational costs. Effects of implicit solvation on molecular geometries and SPEs, and their overall effects on the prediction accuracy of redox potentials are analyzed in view of computational cost versus prediction accuracy, which outlines the best choice of methods corresponding to a desired level of accuracy. The modular computational approach is applicable for accelerating the virtual studies on functional quinones and the respective discovery of candidate compounds for energy storage.

Reduced-cost supercell approach for computing accurate phonon density of states in organic crystals

Phonon contributions to organic crystal structures and thermochemical properties can be significant, but computing a well-converged phonon density of states with lattice dynamics and periodic density functional theory (DFT) is often computationally expensive due to the need for large supercells. Using semi-empirical methods like density functional tight binding (DFTB) instead of DFT can reduce the computational costs dramatically, albeit with noticeable reductions in accuracy. This work proposes approximating the phonon density of states via a relatively inexpensive DFTB supercell treatment of the phonon dispersion that is then corrected by shifting the individual phonon modes according to the difference between the DFT and DFTB phonon frequencies at the Γ-point. The acoustic modes are then computed at the DFT level from the elastic constants. In several small-molecule crystal test cases, this combined approach reproduces DFT thermochemistry with kJ/mol accuracy and 1-2 orders of magnitude less computational effort. Finally, this approach is applied to computing the free energy differences between the five crystal polymorphs of oxalyl dihydrazide.

A quantum-based molecular dynamics study of the ICM-102/HNO3 host-guest reaction at high temperatures

The contradiction between energy and safety of explosives is better balanced by the host-guest inclusion strategy. Understanding the reaction mechanism of the host-guest explosive is necessary. To deeply analyze the role of the small guest molecules in the host-guest system, a quantum-based molecular dynamics method was used to calculate the initial decomposition reaction of the new host-guest explosive ICM-102/HNO3 against the pure ICM-102 at several high temperatures. The incorporation of HNO3 had no significant influence on the initial decomposition step of ICM-102. Conversely, the earliest intramolecular hydrogen transfer reaction is delayed partly because the H and O atoms of HNO3 connect with the O and H atoms of ICM-102, respectively. As the reaction proceeds, guest molecules get heavily involved in the reaction and increase the reaction rate. The generation rate and quantity of the small oxidizing molecules in the final product were increased significantly in the ICM-102/HNO3 system. These mechanisms revealed that HNO3 molecules inhibit the early stages of the initial decomposition of ICM-102 to some extent, and play an important role in accelerating the decomposition subsequently.

Adventures in DFTB: Toward an Automatic Parameterization Scheme

As we push forward on understanding the fate of chemicals in the environment, we need a method that will allow for the simulation of the inherent heterogeneity. Density functional tight binding (DFTB) is a methodology that allows for a detailed electronic description and would be ideal for this problem. While many parameters can be derived directly from DFT, empirical parameters still exist in the confinement and repulsion potentials. In this manuscript, we examine these potentials and present solutions that will minimize the degree of empiricism. Our results show that it is possible to construct confinement potentials from examining the atomic radial wavefunctions. Moreover, we found that the heterogeneous repulsion potentials can be derived from using only homogeneous repulsion curves.

Free-energy landscape of polymer-crystal polymorphism

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