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.

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.

Transferability of Buckingham Parameters for Short-Range Repulsion between Topological Atoms.

The repulsive part of the Buckingham potential, with parameters A and B, can be used to model deformation energies and steric energies. Both are calculated using the interacting quantum atom energy decomposition scheme where the latter is obtained from the former by a charge-transfer-based energy correction. These energies relate to short-range interactions, specifically the deformation of electron density and steric hindrance, respectively, when topological atoms approach each other. In this work, we calculate and fit the energies of carbonyl carbon, carbonyl oxygen, and, where possible, amine nitrogen atoms to the repulsive part of the Buckingham potential for 26 molecules. We find that while the steric energies of all atom pairs studied display exponential behavior with respect to distance, some deformation energy data do not. The obtained parameters are shown to be transferable by calculating root-mean-square errors of fitted potentials with respect to energy data of the same atom in, as far as possible, all other molecules from our data set. We observed that 36% and 10% of these errors were smaller than 4 kJ mol-1 for steric and deformation energy, respectively. Thus, we find that steric energy parameters are more transferable than deformation energy parameters. Finally, we speculate about the physical meaning of the A and B parameters and the implications of the strong exponential and exponential-linear piecewise relationships that we observe between them.

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.

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.

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

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.

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.

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.

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.

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.

On the many-body nature of intramolecular forces in FFLUX and its implications.

FFLUX is a biomolecular force field under construction, based on Quantum Chemical Topology (QCT) and machine learning (kriging), with a minimalistic and physically motivated design. A detailed analysis of the forces within the kriging models as treated in FFLUX is presented, taking as a test example a liquid water model. The energies of topological atoms are modeled as 3Natoms -6 dimensional potential energy surfaces, using atomic local frames to represent the internal degrees of freedom. As a result, the forces within the kriging models in FFLUX are inherently N-body in nature where N refers to Natoms . This provides a fuller picture that is closer to a true quantum mechanical representation of interactions between atoms. The presented computational example quantitatively showcases the non-negligible (as much as 9%) three-body nature of bonded forces and angular forces in a water molecule. We discuss the practical impact on the pressure calculation with N-body forces and periodic boundary conditions (PBC) in molecular dynamics, as opposed to classical force fields with two-body forces. The equivalence between the PBC-related correction terms in the general virial equation is shown mathematically.

Electrostatics Explains the Reverse Lewis Acidity of BH3 and Boron Trihalides: Infrared Intensities and a Relative Energy Gradient (REG) Analysis of IQA Energies.

The reaction path for the formation of BX3-NH3 (X = H, F, Cl, Br) complexes was divided into two processes: (i) rehybridization of the acid while adopting a pyramidal geometry, and (ii) the complex formation from the pyramidal geometries of the acid and base. The interacting quantum atom (IQA) method was used to investigate the Lewis acidity trend of these compounds. This topological analysis suggests that the boron-halogen bond exhibits a considerable degree of ionicity. A relative energy gradient (REG) analysis on IQA energies indicates that the acid-base complex formation is highly dependent on electrostatic energy. With increasing halogen electronegativity, a higher degree of ionicity of the B-X is observed, causing an increase in the absolute value of X and B charges. This increases not only the attractive electrostatic energy between the acid and base but also enhances the repulsive energy. The latter is the main factor behind the acidity trend exhibited by trihalides. Changes in geometry are relevant only for complexes where BH3 acts as an acid, where lower steric hindrance facilitates the adoption of the pyramidal geometry observed in the complex. The CCTDP analysis shows that infrared intensities of BX3-NH3 are determined mostly by the atomic charges and not by the charge transfer or polarization. The opposite is observed in covalent analogues.

Enhancing Carbon Acid pKa Prediction by Augmentation of Sparse Experimental Datasets with Accurate AIBL (QM) Derived Values.

The prediction of the aqueous pKa of carbon acids by Quantitative Structure Property Relationship or cheminformatics-based methods is a rather arduous problem. Primarily, there are insufficient high-quality experimental data points measured in homogeneous conditions to allow for a good global model to be generated. In our computationally efficient pKa prediction method, we generate an atom-type feature vector, called a distance spectrum, from the assigned ionisation atom, and learn coefficients for those atom-types that show the impact each atom-type has on the pKa of the ionisable centre. In the current work, we augment our dataset with pKa values from a series of high performing local models derived from the Ab Initio Bond Lengths method (AIBL). We find that, in distilling the knowledge available from multiple models into one general model, the prediction error for an external test set is reduced compared to that using literature experimental data alone.

A FFLUX Water Model: Flexible, Polarizable and with a Multipolar Description of Electrostatics.

Key to progress in molecular simulation is the development of advanced models that go beyond the limitations of traditional force fields that employ a fixed, point charge-based description of electrostatics. Taking water as an example system, the FFLUX framework is shown capable of producing models that are flexible, polarizable and have a multipolar description of the electrostatics. The kriging machine-learning methods used in FFLUX are able to reproduce the intramolecular potential energy surface and multipole moments of a single water molecule with chemical accuracy using as few as 50 training configurations. Molecular dynamics simulations of water clusters (25-216 molecules) using the new FFLUX model reveal that incorporating charge-quadrupole, dipole-dipole, and quadrupole-charge interactions into the description of the electrostatics results in significant changes to the intermolecular structuring of the water molecules. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Creating Gaussian process regression models for molecular simulations using adaptive sampling.

FFLUX is a new force field that combines the accuracy of quantum mechanics with the speed of force fields, without any link to the architecture of classical force fields. This force field is atom-focused and adopts the parameter-free topological atom from Quantum Chemical Topology (QCT). FFLUX uses Gaussian process regression (also known as kriging) models to make predictions of atomic properties, which in this work are atomic energies according to QCT's interacting quantum atom approach. Here, we report the adaptive sampling technique maximum expected prediction error to create data-compact, efficient, and accurate kriging models (sub-kJ mol-1 for water, ammonia, methane, and methanol and sub-kcal mol-1 for N-methylacetamide). The models cope with large molecular distortions and are ready for use in molecular simulation. A brand new press-one-button Python pipeline, called ICHOR, carries out the training.