Determining the Factors accounting for Reaction Selectivity A Relative Energy Gradient - Interacting Quantum Atoms and Natural Bonding Orbitals Study.

Identifying the main physicochemical properties accounting for the course of a reaction is of utmost importance to rationalize chemical syntheses. To this aim, the relative energy gradient (REG) method is an appealing approach because it is an unbiased and automatic process to extract the most relevant pieces of energy information. Initially formulated within the interacting quantum atoms (IQA) framework for a single reaction, here we extend the REG method to natural bond orbitals (NBO) analysis and to the case of two competitive processes. This development enables the determination of the driving forces of any chemical selectivity. We illustrate the extended REG method on the case study of ring opening in cyclobutenes, which is an important instance of the so-called torquoselectivity.

An Unsupervised Machine Learning Approach for the Automatic Construction of Local Chemical Descriptors.

Condensing the many physical variables defining a chemical system into a fixed-size array poses a significant challenge in the development of chemical Machine Learning (ML). Atom Centered Symmetry Functions (ACSFs) offer an intuitive featurization approach by means of a tedious and labor-intensive selection of tunable parameters. In this work, we implement an unsupervised ML strategy relying on a Gaussian Mixture Model (GMM) to automatically optimize the ACSF parameters. GMMs effortlessly decompose the vastness of the chemical and conformational spaces into well-defined radial and angular clusters, which are then used to build tailor-made ACSFs. The unsupervised exploration of the space has demonstrated general applicability across a diverse range of systems, spanning from various unimolecular landscapes to heterogeneous databases. The impact of the sampling technique and temperature on space exploration is also addressed, highlighting the particularly advantageous role of high-temperature Molecular Dynamics (MD) simulations. The reliability of the resulting features is assessed through the estimation of the atomic charges of a prototypical capped amino acid and a heterogeneous collection of CHON molecules. The automatically constructed ACSFs serve as high-quality descriptors, consistently yielding typical prediction errors below 0.010 electrons bound for the reported atomic charges. Altering the spatial distribution of the functions with respect to the cluster highlights the critical role of symmetry rupture in achieving significantly improved features. More specifically, using two separate functions to describe the lower and upper tails of the cluster results in the best performing models with errors as low as 0.006 electrons. Finally, the effectiveness of finely tuned features was checked across different architectures, unveiling the superior performance of Gaussian Process (GP) models over Feed Forward Neural Networks (FFNNs), particularly in low-data regimes, with nearly a 2-fold increase in prediction quality. Altogether, this approach paves the way toward an easier construction of local chemical descriptors, while providing valuable insights into how radial and angular spaces should be mapped. Finally, this work opens the possibility of encoding many-body information beyond angular terms into upcoming ML features.

Probing Non-Covalent Interactions through Molecular Balances: A REG-IQA Study.

The interaction energies of two series of molecular balances (1-X with X = H, Me, OMe, NMe2 and 2-Y with Y = H, CN, NO2, OMe, NMe2) designed to probe carbonyl…carbonyl interactions were analysed at the B3LYP/6-311++G(d,p)-D3 level of theory using the energy partitioning method of Interacting Quantum Atoms/Fragments (IQA/IQF). The partitioned energies are analysed by the Relative Energy Gradient (REG) method, which calculates the correlation between these energies and the total energy of a system, thereby explaining the role atoms have in the energetic behaviour of the total system. The traditional "back-of-the-envelope" open and closed conformations of molecular balances do not correspond to those of the lowest energy. Hence, more care needs to be taken when considering which geometries to use for comparison with the experiment. The REG-IQA method shows that the 1-H and 1-OMe balances behave differently to the 1-Me and 1-NMe2 balances because the latter show more prominent electrostatics between carbonyl groups and undergoes a larger dihedral rotation due to the bulkiness of the functional groups. For the 2-Y balance, REG-IQA shows the same behaviour across the series as the 1-H and 1-OMe balances. From an atomistic point of view, the formation of the closed conformer is favoured by polarisation and charge-transfer effects on the amide bond across all balances and is counterbalanced by a de-pyramidalisation of the amide nitrogen. Moreover, focusing on the oxygen of the amide carbonyl and the α-carbon of the remaining carbonyl group, electrostatics have a major role in the formation of the closed conformer, which goes against the well-known n-π* interaction orbital overlap concept. However, REG-IQF shows that exchange-correlation energies overtake electrostatics for all the 2-Y balances when working with fragments around the carbonyl groups, while they act on par with electrostatics for the 1-OMe and 1-NMe2. REG-IQF also shows that exchange-correlation energies in the 2-Y balance are correlated to the inductive electron-donating and -withdrawing trends on aromatic groups. We demonstrate that methods such as REG-IQA/IQF can help with the fine-tuning of molecular balances prior to the experiment and that the energies that govern the probed interactions are highly dependent on the atoms and functional groups involved.

FFLUX molecular simulations driven by atomic Gaussian process regression models.

Machine learning (ML) force fields are revolutionizing molecular dynamics (MD) simulations as they bypass the computational cost associated with ab initio methods but do not sacrifice accuracy in the process. In this work, the GPyTorch library is used to create Gaussian process regression (GPR) models that are interfaced with the next-generation ML force field FFLUX. These models predict atomic properties of different molecular configurations that appear in a progressing MD simulation. An improved kernel function is utilized to correctly capture the periodicity of the input descriptors. The first FFLUX molecular simulations of ammonia, methanol, and malondialdehyde with the updated kernel are performed. Geometry optimizations with the GPR models result in highly accurate final structures with a maximum root-mean-squared deviation of 0.064 Å and sub-kJ mol-1 total energy predictions. Additionally, the models are tested in 298 K MD simulations with FFLUX to benchmark for robustness. The resulting energy and force predictions throughout the simulation are in excellent agreement with ab initio data for ammonia and methanol but decrease in quality for malondialdehyde due to the increased system complexity. GPR model improvements are discussed, which will ensure the future scalability to larger systems.

Application of the FFLUX Force Field to Molecular Crystals: A Study of Formamide.

In this work, we present the first application of the quantum chemical topology force field FFLUX to the solid state. FFLUX utilizes Gaussian process regression machine learning models trained on data from the interacting quantum atom partitioning scheme to predict atomic energies and flexible multipole moments that change with geometry. Here, the ambient (α) and high-pressure (ß) polymorphs of formamide are used as test systems and optimized using FFLUX. Optimizing the structures with increasing multipolar ranks indicates that the lattice parameters of the α phase differ by less than 5% to the experimental structure when multipole moments up to the quadrupole are used. These differences are found to be in line with the dispersion-corrected density functional theory. Lattice dynamics calculations are also found to be possible using FFLUX, yielding harmonic phonon spectra comparable to dispersion-corrected DFT while enabling larger supercells to be considered than is typically possible with first-principles calculations. These promising results indicate that FFLUX can be used to accurately determine properties of molecular solids that are difficult to access using DFT, including the structural dynamics, free energies, and properties at finite temperature.

Interacting Quantum Atoms and Multipolar Electrostatic Study of XH···π Interactions.

The interaction energies of nine XH···π (X = C, N, and O) benzene-containing van der Waals complexes were analyzed, at the atomic and fragment levels, using QTAIM multipolar electrostatics and the energy partitioning method interacting quantum atoms/fragment (IQA/IQF). These descriptors were paired with the relative energy gradient method, which solidifies the connection between quantum mechanical properties and chemical interpretation. This combination provides a precise understanding, both qualitative and quantitative, of the nature of these interactions, which are ubiquitous in biochemical systems. The formation of the OH···π and NH···π systems is electrostatically driven, with the Qzz component of the quadrupole moment of the benzene carbons interacting with the charges of X and H in XH. There is the unexpectedly intramonomeric role of X-H (X = O, N) where its electrostatic energy helps the formation of the complex and its covalent energy thwarts it. However, the CH···π interaction is governed by exchange-correlation energies, thereby establishing a covalent character, as opposed to the literature's designation as a noncovalent interaction. Moreover, dispersion energy is relevant, statically and in absolute terms, but less relevant compared to other energy components in terms of the formation of the complex. Multipolar electrostatics are similar across all systems.

Aza-Michael Addition in Explicit Solvent: A Relative Energy Gradient-Interacting Quantum Atoms Study.

Aza-Michael additions are key reactions in organic synthesis. We investigate, from a theoretical and computational point of view, several examples ranging from weak to strong electrophiles in dimethylsulfoxide treated as explicit solvent. We use the REG-IQA method, which is a quantum topological energy decomposition (Interacting Quantum Atoms, IQA) coupled to a chemical-interpretation calculator (Relative Energy Gradient, REG). We focus on the rate-limiting addition step in order to unravel the different events taking place in this step, and understand the influence of solvent on the reaction, with an eye on predicting the Mayr electrophilicity. For the first time, a link is established between an REG-IQA analysis and experimental values.

How to Compute Atomistic Insight in DFT Clusters: The REG-IQA Approach.

The relative energy gradient (REG) method is paired with the topological energy partitioning method interacting quantum atoms (IQA), as REG-IQA, to provide detailed and unbiased knowledge on the intra- and interatomic interactions. REG operates on a sequence of geometries representing a dynamical change of a system. Its recent application to peptide hydrolysis of the human immunodeficiency virus-1 (HIV-1) protease (PDB code: 4HVP) has demonstrated its full potential in recovering reaction mechanisms and through-space electrostatic and exchange-correlation effects, making it a compelling tool for analyzing enzymatic reactions. In this study, the computational efficiency of the REG-IQA method for the 133-atom HIV-1 protease quantum mechanical system is analyzed in every detail and substantially improved by means of three different approaches. The first approach of smaller integration grids for IQA integrations reduces the computational overhead by about a factor of 3. The second approach uses the line-simplification Ramer-Douglas-Peucker (RDP) algorithm, which outputs the minimal number of geometries necessary for the REG-IQA analysis for a predetermined root mean squared error (RMSE) tolerance. This cuts the computational time of the whole REG analysis by a factor of 2 if an RMSE of 0.5 kJ/mol is considered. The third approach consists of a "biased" or "unbiased" selection of a specific subset of atoms of the whole initial quantum mechanical model wave-function, which results in more than a 10-fold speed-up per geometry for the IQA calculation, without deterioration of the outcome of the REG-IQA analysis. Finally, to show the capability of these approaches, the findings gathered from the HIV-1 protease system are also applied to a different system named haloalcohol dehalogenase (HheC). In summary, this study takes the REG-IQA method to a computationally feasible and highly accurate level, making it viable for the analysis of a multitude of enzymatic systems.

A combined BET and IQA-REG study of the activation energy of non-polar zw-type [3+2] cycloaddition reactions.

A combined Bonding Evolution Theory (BET) and Interacting Quantum Atoms-Relative Energy Gradient (IQA-REG) study is carried out on a non-polar zw-type [3+2] cycloaddition (32CA) reaction. BET is the joint use of Catastrophe Theory and the topology of the Electron Localization Function (ELF) to characterise molecular mechanisms, while IQA is a quantum topological energy partitioning method and REG is a method to compute chemical insight at atomistic level, usually in connection with energy. This 32CA reaction involves the simplest nitrone with ethylene and has been studied here at B3LYP/6-311G(d,p) level within the context of Molecular Electron Density Theory (MEDT), which is based on the idea that changes in electron density, and not molecular orbital interactions, are responsible for chemical reactivity. We aim to determine the origin of the high activation energy of 32CA reactions involving zwitterionic three-atom-components. The BET study and IQA-REG method are applied to the overall activation energy path. While BET suggests that the barrier is mainly associated with the rupture of the nitrone CN double bond, IQA-REG indicates that it is mainly related to the rupture of the ethylene CC double bond. The present study shows that activation energies can be accurately and easily described by IQA-REG, and its complementary use with BET helps achieving a more detailed description of molecular mechanisms.

Construction of a Gaussian Process Regression Model of Formamide for Use in Molecular Simulations.

FFLUX, a novel force field based on quantum chemical topology, can perform molecular dynamics simulations with flexible multipole moments that change with geometry. This is enabled by Gaussian process regression machine learning models, which accurately predict atomic energies and multipole moments up to the hexadecapole. We have constructed a model of the formamide monomer at the B3LYP/aug-cc-pVTZ level of theory capable of sub-kJ mol-1 accuracy, with the maximum prediction error for the molecule being 0.8 kJ mol-1. This model was used in FFLUX simulations along with Lennard-Jones parameters to successfully optimize the geometry of formamide dimers with errors smaller than 0.1 Šcompared to those obtained with D3-corrected B3LYP/aug-cc-pVTZ. Comparisons were also made to a force field constructed with static multipole moments and Lennard-Jones parameters. FFLUX recovers the expected energy ranking of dimers compared to the literature, and changes in C═O and C-N bond lengths associated with hydrogen bonding were found to be consistent with density functional theory.

Gaussian Process Regression Models for Predicting Atomic Energies and Multipole Moments.

Developing a force field is a difficult task because its design is typically pulled in opposite directions by speed and accuracy. FFLUX breaks this trend by utilizing Gaussian process regression (GPR) to predict, at ab initio accuracy, atomic energies and multipole moments as obtained from the quantum theory of atoms in molecules (QTAIM). This work demonstrates that the in-house FFLUX training pipeline can generate successful GPR models for six representative molecules: peptide-capped glycine and alanine, glucose, paracetamol, aspirin, and ibuprofen. The molecules were sufficiently distorted to represent configurations from an AMBER-GAFF2 molecular dynamics run. All internal degrees of freedom were covered corresponding to 93 dimensions in the case of the largest molecule ibuprofen (33 atoms). Benefiting from active learning, the GPR models contain only about 2000 training points and return largely sub-kcal mol-1 prediction errors for the validation sets. A proof of concept has been reached for transferring the model produced through active learning on one atomic property to that of the remaining atomic properties. The prediction of electrostatic interaction can be assessed at the intermolecular level, and the vast majority of interactions have a root-mean-square error of less than 0.1 kJ mol-1 with a maximum value of ∼1 kJ mol-1 for a glycine and paracetamol dimer.

Linking the Interatomic Exchange-Correlation Energy to Experimental J-Coupling Constants.

The main aim of the current work is to find an experimental connection to the interatomic exchange-correlation energy as defined by the energy decomposition method Interacting Quantum Atoms (IQA). A suitable candidate as (essentially) experimental quantity is the nuclear magnetic resonance (NMR) J-coupling constant denoted 3J(H,H'), which a number of previous studies showed to correlate well with QTAIM's delocalization index (DI), which is essentially a bond order. Inspired by Karplus equations, here, we investigate correlations between 3J(H,H') and a relevant dihedral angle in six simple initial compounds of the shape H3C-YHn (Y = C, N, O, Si, P, and S), N-methylacetamide (as prototype of the peptide bond), and five peptide-capped amino acids (Gly, Ala, Val, Ile, and Leu) because of the protein direction of the force field FFLUX. In conclusion, except for methanol, the inter-hydrogen exchange-correlation energy Vxc(H,H') makes the best contact with experiment, through 3J(H,H'), when multiplied with the internuclear distance RHH'.

The IQA Energy Partition in a Drug Design Setting: A Hepatitis C Virus RNA-Dependent RNA Polymerase (NS5B) Case Study.

The interaction of the thumb site II of the NS5B protein of hepatitis C virus and a pair of drug candidates was studied using a topological energy decomposition method called interacting quantum atoms (IQA). The atomic energies were then processed by the relative energy gradient (REG) method, which extracts chemical insight by computation based on minimal assumptions. REG reveals the most important IQA energy contributions, by atom and energy type (electrostatics, sterics, and exchange-correlation), that are responsible for the behaviour of the whole system, systematically from a short-range ligand-pocket interaction until a distance of approximately 22 Å. The degree of covalency in various key interatomic interactions can be quantified. No exchange-correlation contribution is responsible for the changes in the energy profile of both pocket-ligand systems investigated in the ligand-pocket distances equal to or greater than that of the global minimum. Regarding the hydrogen bonds in the system, a "neighbour effect" was observed thanks to the REG method, which states that a carbon atom would rather not have its covalent neighbour oxygen form a hydrogen bond. The combination of IQA and REG enables the automatic identification of the pharmacophore in the ligands. The coarser Interacting Quantum Fragments (IQF) enables the determination of which amino acids of the pocket contribute most to the binding and the type of energy of said binding. This work is an example of the contribution topological energy decomposition methods can make to fragment-based drug design.

Producing chemically accurate atomic Gaussian process regression models by active learning for molecular simulation.

Machine learning is becoming increasingly more important in the field of force field development. Never has it been more vital to have chemically accurate machine learning potentials because force fields become more sophisticated and their applications expand. In this study a method for developing chemically accurate Gaussian process regression models is demonstrated for an increasingly complex set of molecules. This work is an extension to previous work showing the progression of the active learning technique in producing more accurate models in much less CPU time than ever before. The per-atom active learning approach has unlocked the potential to generate chemically accurate models for molecules such as peptide-capped glycine.

Binding Energy Partition of Promising IRAK-4 Inhibitor (Zimlovisertib) for the Treatment of COVID-19 Pneumonia.

The technique of Fragment-Based Drug Design (FBDD) considers the interactions of different moieties of molecules with biological targets for the rational construction of potential drugs. One basic assumption of FBDD is that the different functional groups of a ligand interact with a biological target in an approximately additive, that is, independent manner. We investigated the interactions of different fragments of ligands and Interleukin-1 Receptor-Associated Kinase 4 (IRAK-4) throughout the FBDD design of Zimlovisertib, a promising anti-inflammatory, currently in trials to be used for the treatment of COVID-19 pneumonia. We utilised state-of-the-art methods of wave function analyses mainly the Interacting Quantum Atoms (IQA) energy partition for this purpose. By means of IQA, we assessed the suitability of every change to the ligand in the five stages of FBDD which led to Zimlovisertib on a quantitative basis. We determined the energetics of the interaction of different functional groups in the ligands with the IRAK-4 protein target and thereby demonstrated the adequacy (or lack thereof) of the changes made across the design of this drug. This analysis permits to verify whether a given alteration of a prospective drug leads to the intended tuning of non-covalent interactions with its protein objective. Overall, we expect that the methods exploited in this paper will prove valuable in the understanding and control of chemical modifications across FBDD processes.

An Interacting Quantum Atoms (IQA) and Relative Energy Gradient (REG) Analysis of the Anomeric Effect.

The explanation of the anomeric effect in terms of underlying quantum properties is still controversial almost 70 years after its introduction. Here, we use a method called Relative Energy Gradient (REG), which is able to compute chemical insight with a view to explaining the anomeric effect. REG operates on atomic energy contributions generated by the quantum topological energy decomposition Interacting Quantum Atoms (IQA). Based on the case studies of dimethoxymethane and 2-fluorotetrahydropyran, we show that the anomeric effect is electrostatic in nature rather than governed by hyperconjugation.

Application of Quantum Chemical Topology Force Field FFLUX to Condensed Matter Simulations: Liquid Water.

We present here the first application of the quantum chemical topology force field FFLUX to condensed matter simulations. FFLUX offers many-body potential energy surfaces learnt exclusively from ab initio data using Gaussian process regression. FFLUX also includes high-rank, polarizable multipole moments (up to quadrupole moments in this work) that are learnt from the same ab initio calculations as the potential energy surfaces. Many-body effects (where a body is an atom) and polarization are captured by the machine learning models. The choice to use machine learning in this way allows the force field's representation of reality to be improved (e.g., by including higher order many-body effects) with little detriment to the computational scaling of the code. In this manner, FFLUX is inherently future-proof. The "plug and play" nature of the machine learning models also ensures that FFLUX can be applied to any system of interest, not just liquid water. In this work we study liquid water across a range of temperatures and compare the predicted bulk properties to experiment as well as other state-of-the-art force fields AMOEBA(+CF), HIPPO, MB-Pol and SIBFA21. We find that FFLUX finds a place amongst these.

Non-covalent interactions from a Quantum Chemical Topology perspective.

About half a century after its little-known beginnings, the quantum topological approach called QTAIM has grown into a widespread, but still not mainstream, methodology of interpretational quantum chemistry. Although often confused in textbooks with yet another population analysis, be it perhaps an elegant but somewhat esoteric one, QTAIM has been enriched with about a dozen other research areas sharing its main mathematical language, such as Interacting Quantum Atoms (IQA) or Electron Localisation Function (ELF), to form an overarching approach called Quantum Chemical Topology (QCT). Instead of reviewing the latter's role in understanding non-covalent interactions, we propose a number of ideas emerging from the full consequences of the space-filling nature of topological atoms, and discuss how they (will) impact on interatomic interactions, including non-covalent ones. The architecture of a force field called FFLUX, which is based on these ideas, is outlined. A new method called Relative Energy Gradient (REG) is put forward, which is able, by computation, to detect which fragments of a given molecular assembly govern the energetic behaviour of this whole assembly. This method can offer insight into the typical balance of competing atomic energies both in covalent and non-covalent case studies. A brief discussion on so-called bond critical points is given, highlighting concerns about their meaning, mainly in the arena of non-covalent interactions.

Flexible multipole moments in smooth particle mesh Ewald.

The smooth particle mesh Ewald sum is extended with additional force terms that arise from the so-called flexible multipole moments. These are multipole moments (of any rank) that depend explicitly on atomic positions in some local environment that can be made arbitrarily large. By introducing explicit dependence on atomic positions, flexible multipole moments are polarized by their local environment, allowing both intramolecular and intermolecular polarizations to be captured. Multipolar torques are discussed in detail, and it is shown that they arise naturally in the presented framework. Furthermore, we give details of how we validated our implementation of the flexible smooth particle mesh Ewald sum by considering two mathematical limits of the smooth particle mesh Ewald summation.

Towards an atomistic understanding of polymorphism in molecular solids.

Understanding and controlling polymorphism in molecular solids is a major unsolved problem in crystal engineering. While the ability to calculate accurate lattice energies with atomistic modelling provides valuable insight into the associated energy scales, existing methods cannot connect energy differences to the delicate balances of intra- and intermolecular forces that ultimately determine polymorph stability ordering. We report herein a protocol for applying Quantum Chemical Topology (QCT) to study the key intra- and intermolecular interactions in molecular solids, which we use to compare the three known polymorphs of succinic acid including the recently-discovered γ form. QCT provides a rigorous partitioning of the total energy into contributions associated with topological atoms, and a quantitative and chemically intuitive description of the intra- and intermolecular interactions. The newly-proposed Relative Energy Gradient (REG) method ranks atomistic energy terms (steric, electrostatic and exchange) by their importance in constructing the total energy profile for a chemical process. We find that the conformation of the succinic acid molecule is governed by a balance of large and opposing electrostatic interactions, while the H-bond dimerisation is governed by a combination of electrostatics and sterics. In the solids, an atomistic energy balance emerges that governs the contraction, towards the equilibrium geometry, of a molecular cluster representing the bulk crystal. The protocol we put forward is as general as the capabilities of the underlying quantum-mechanical model and it can provide novel perspectives on polymorphism in a wide range of chemical systems.