Automated Discovery and Refinement of Reactive Molecular Dynamics Pathways

Facilitating Transition State Search with Minimal Conformational Sampling Using Reaction Graph

Elucidating transition states (TSs) is crucial for understanding chemical reactions. The reliability of traditional TS search approaches depends on input conformations that require significant effort to prepare. Previous automated methods for generating input reaction conformations typically involve extensive exploration of a large conformational space. Such exhaustive search can be complicated by the rapid growth of the conformational space, especially for reactions involving many rotatable bonds, multiple reacting molecules, and numerous bond formations and dissociations. To address this problem, we propose a new approach that generates reaction conformations for TS searches with minimal reliance on sampling. This method constructs a pseudo-TS structure based on a reaction graph containing bond formation and dissociation information and modifies it to produce reactant and product conformations. Tested on three different benchmarks, our method consistently generated suitable conformations without necessitating extensive sampling, demonstrating its potential to significantly improve the applicability of automated TS searches. This approach offers a valuable tool for a broad range of applications such as reaction mechanism analysis and network exploration.

Chemical Reaction Networks from Scratch with Reaction Prediction and Kinetics-Guided Exploration

Algorithmic reaction explorations based on transition state searches can now routinely predict relatively short reaction sequences involving small molecules. However, applying these algorithms to deeper chemical reaction network (CRN) exploration still requires the development of more efficient and accurate exploration policies. Here, an exploration algorithm, which we name yet another kinetic strategy (YAKS), is demonstrated that uses microkinetic simulations of the nascent network to achieve cost-effective, deep network exploration. Key features of the algorithm are the automatic incorporation of bimolecular reactions between network intermediates, compatibility with short-lived but kinetically important species, and incorporation of rate uncertainty into the exploration policy. In validation case studies of glucose pyrolysis, the algorithm rediscovers reaction pathways previously discovered by heuristic exploration policies and elucidates new reaction pathways for experimentally obtained products. The resulting CRN is the first to connect all major experimental pyrolysis products to glucose. Additional case studies are presented that investigate the role of reaction rules, rate uncertainty, and bimolecular reactions. These case studies show that naïve exponential growth estimates can vastly overestimate the actual number of kinetically relevant pathways in the physical reaction networks. In light of this, further improvements in exploration policies and reaction prediction algorithms make it feasible that CRNs might soon be routinely predictable in some contexts.

Automated search of minimum-energy conical intersections with projected metadynamics

We present a new method for the automated search of minimum-energy conical intersections (MECIs) based on metadynamics. In this method, two independent forces are constructed and projected into the minimization subspace and the constraint subspace, respectively. One force is directed toward the minimum-energy point, while the other is directed toward the conical intersection seam. The root-mean-square deviation based bias potential is added to the potential energy surface to force the structure escape from the already explored regions. The additional constraint function is used to enable the structure reach different intersection seams. This method can be used for systematically and automatically searching MECIs or exploring conical intersection seams. Compared to the penalty function-based metadynamics method, this new method is more effective and stable in searching MECIs. Furthermore, this method can be combined with any kind of constraint, whether geometric or non-geometric, making it a generalized tool for the automated search of constrained minimum.

Acceleration of Diffusion in Ab Initio Nanoreactor Molecular Dynamics and Application to Hydrogen Sulfide Oxidation

The computational description of chemical reactivity can become extremely complex when multiple different reaction products and intermediates come into play, forming a chemical reaction network. Therefore, computational methods for the automated construction of chemical reaction networks have been developed in the last decades. One of these methods, ab initio nanoreactor molecular dynamics (NMD), is based on external forces enhancing reactivity by e.g., periodically compressing the system and allowing it to relax. However, during the relaxation process, a significant simulation time is required to allow energy to dissipate and molecules to diffuse, making this part of the NMD simulation computationally intensive. This work aims to improve NMD by accelerating the diffusion process in the relaxation phase. We systematically investigate the speedup of reaction discovery gained by diffusion acceleration, leading to a factor of up to 28 in discovery frequency. Diffusion-accelerated nanoreactor molecular dynamics (DA-NMD) is then used to construct a reaction network of hydrogen sulfide oxidation under atmospheric conditions, where reactions are automatically detected by a change in the bond order and bond distance. A reaction network of 108 molecular species and 399 elementary reactions was constructed starting from hydrogen sulfide, hydroxy radicals, and molecular oxygen covering a broad variety of sulfur-oxygen chemistry and oxidation states of the sulfur atom ranging from -II to +VI.

Ab Initio Molecular Dynamics Study of Trivalent Rare Earth Rich Borate Glasses: Structural Insights and Formation Mechanisms

In this work, trivalent rare earth (RE)-rich borate glasses (30 La2O3-70 B2O3, 50 La2O3-50 B2O3, 60 La2O3-40 B2O3, and 50 Y2O3-50 B2O3) were modeled using ab initio molecular dynamics (AIMD) simulations through the melt-quenching route. It was found that the AIMD-derived structures reproduced the experimental structure factors and 11B solid-state nuclear magnetic resonance data. Isolated borate units (monomers, dimers, and trimers) terminated with nonbridging oxygen were found in the structures. Polymer units containing four or more boron atoms were identified with and without three-membered boron rings (3-rings). Increasing the proportion of La2O3 in La2O3-B2O3 glasses resulted in an increased number of isolated units, indicating that La3+ acts as a network modifier, breaking the borate glass network. The formation of these units via the melt-quenching process was detected by labeling boron species at each AIMD step from 1500 to 300 K. Representation with transition matrices clarified the specific reaction routes, leading to the formation of isolated boron units in solid glass. A key finding is the stabilization of polymer units involving 3-ring formation. The formation of isolated units is achieved through the reaction of polymers without 3-rings. The RE coordination structure was thoroughly analyzed from the perspective of shape and symmetry. Reference structures derived from the solution of the Thomson problem were compared to the AIMD-derived coordination structures and crystalline LaBO3 and YBO3. The results highlight the specificity of the Y coordination structure with 3-rings in YBO3, which is not observed in RE borate glasses. The analytical approaches and interpretations used in this study provide insights into the diverse coordination structures of glasses containing heavy elements other than REs.

Predicting Collision-Induced-Dissociation Tandem Mass Spectra (CID-MS/MS) Using Ab Initio Molecular Dynamics

Compound identification is at the center of metabolomics, usually by comparing experimental mass spectra against library spectra. However, most compounds are not commercially available to generate library spectra. Hence, for such compounds, MS/MS spectra need to be predicted. Machine learning and heuristic models have largely failed except for lipids. Here, quantum chemistry software can be used to predict mass spectra. However, quantum chemistry predictions for collision induced dissociation (CID) mass spectra in LC-MS/MS are rare. We present the CIDMD (Collision-Induced Dissociation via Molecular Dynamics) framework to model CID-based MS/MS spectra. It uses first-principles molecular dynamics (MD) to simulate the physical process of molecular collisions in CID tandem mass spectrometry. First, molecular ions are constructed at specific protonation sites. Using density functional theory, these protonated ions are targeted by argon collider gas atoms at user-specified velocities. Subsequent bond breakages are simulated over time for at least 1,000 fs. Each simulation is repeated multiple times from various collisional directions. Fragmentations are accumulated over those repeated collisions to generate CIDMD in silico mass spectra. Twelve small metabolites (<205 Da) were selected to test the accuracy of this framework in comparison to experimental MS/MS spectra. When testing different protomers, collider velocities, number of simulations, simulation time and impact factor b cutoffs, we yielded 261 predicted mass spectra. These in silico spectra resulted in entropy similarity scores of an average 624 ± 189 for all 261 spectra compared to their corresponding experimental spectra, which improved to 828 ± 77 when using optimal parameters of the most probable protomers for 12 molecules. With increasing molecular mass, higher velocities achieved better results. Similarly, different protomers showed large differences in fragmentation; hence, with increasing numbers of protomers and tautomers, the average CIDMD prediction accuracy decreased. Mechanistic details showed that specific fragment ions can be produced from different protomers via multiple fragmentation pathways. We propose that CIDMD is a suitable tool to predict mass spectra of small metabolites like produced by the gut microbiome.

Two P, Ten P, White P, Red P: Mechanistic Exploration of the Oligomerization of Red Phosphorus from Diphosphorus with the Ab Initio Nanoreactor

Phosphorus is critical to humans on many fronts, yet we do not have a mechanistic understanding of some of its most basic transformations and reactions─namely the oligomerization of white phosphorus to red. With heat or under ultraviolet (UV) exposure, it has been experimentally demonstrated that white phosphorus dissociates into diphosphorus units which readily form red phosphorus. However, the mechanism of this process is unknown. The ab initio nanoreactor approach was used to explore the potential energy surface of phosphorus clusters. Density functional theory and metadynamics simulations were used to characterize potential reaction pathways. A mechanism for oligomerization is proposed to take place via diphosphorus additions at π-bonds and weak σ-bonds through three membered ring intermediates. Downhill paths through P6 and P8 clusters eventually result in P10 clusters that can oligomerize into red phosphorus chains. The initial, rate limiting step for this process has an energy barrier of 24.2 kcal/mol.

Reactant Discovery with an Ab Initio Nanoreactor: Exploration of Astrophysical N-Heterocycle Precursors and Formation Pathways

The incorporation of nitrogen atoms into cyclic compounds is essential for terrestrial life; nitrogen-containing (N-)heterocycles make up DNA and RNA nucleobases, several amino acids, B vitamins, porphyrins, and other components of biomolecules. The discovery of these molecules on meteorites with non-terrestrial isotopic abundances supports the hypothesis of exogenous delivery of prebiotic material to early Earth; however, there has been no detection of these species in interstellar environments, indicating that there is a need for greater knowledge of their astrochemical formation and destruction pathways. Here, we present results of simulations of gas-phase pyrrole and pyridine formation from an ab initio nanoreactor, a first-principles molecular dynamics simulation method that accelerates reaction discovery by applying non-equilibrium forces that are agnostic to individual reaction coordinates. Using the nanoreactor in a retrosynthetic mode, starting with the N-heterocycle of interest and a radical leaving group, then considering the discovered reaction pathways in reverse, a rich landscape of N-heterocycle-forming reactivity can be found. Several of these reaction pathways, when mapped to their corresponding minimum energy paths, correspond to novel barrierless formation pathways for pyridine and pyrrole, starting from both detected and hypothesized astrochemical precursors. This study demonstrates how first-principles reaction discovery can build mechanistic knowledge in astrochemical environments as well as in early Earth models such as Titan's atmosphere where N-heterocycles have been tentatively detected.

Reaction Networks Resemble Low-Dimensional Regular Lattices

The computational exploration, manipulation, and design of complex chemical reactions face fundamental challenges related to the high-dimensional nature of potential energy surfaces (PESs) that govern reactivity. Accurately modeling complex reactions is crucial for understanding the chemical processes involved in, for example, organocatalysis, autocatalytic cycles, and one-pot molecular assembly. Our prior research demonstrated that discretizing PESs using heuristics based on bond breaking and bond formation produces a reaction network representation with a low-dimensional structure (metric space). We now find that these stoichiometry-preserving reaction networks possess additional, though approximate, structure and resemble low-dimensional regular lattices with a small amount of random edge rewiring. The heuristics-based discretization thus generates a nonlinear dimensionality reduction by a factor of 10 with an a posteriori error measure (probability of random rewiring). The structure becomes evident through a comparative analysis of CHNO reaction networks of varying stoichiometries against a panel of size-matched generative network models, taking into account their local, metric, and global properties. The generative models include random networks (Erdős-Rényi and bipartite random networks), regular lattices (periodic and nonperiodic), and network models with a tunable level of "randomness" (Watts-Strogatz graphs and regular lattices with random rewiring). The CHNO networks are simultaneously closely matched in all these properties by 3-4-dimensional regular lattices with 10% or less of edges randomly rewired. The effective dimensionality reduction is found to be independent of the system size, stoichiometry, and ruleset, suggesting that search and sampling algorithms for PESs of complex chemical reactions can be effectively leveraged.

A Polarizable QM/MM Model That Combines the State-Averaged CASSCF and AMOEBA Force Field for Photoreactions in Proteins

This study presents the polarizable quantum mechanics/molecular mechanics (QM/MM) embedding of the state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field for the purpose of studying photoreactions in protein environments. We describe two extensions of our previous work that combine SA-CASSCF with AMOEBA water models, allowing it to be generalized to AMOEBA models for proteins and other macromolecules. First, we discuss how our QM/MM model accounts for the discrepancy between the direct and polarization electric fields that arises in the AMOEBA description of intramolecular polarization. A second improvement is the incorporation of link atom schemes to treat instances in which the QM/MM boundary goes through covalent bonds. A single-link atom scheme and double-link atom scheme are considered in this work, and we will discuss how electrostatic interaction, van der Waals interaction, and various kinds of valence terms are treated across the boundary. To test the accuracy of the link atom scheme, we will compare QM/MM with full QM calculations and study how the errors in ground state properties, excited state properties, and excitation energies change when tuning the parameters in the link atom scheme. We will also test the new SA-CASSCF/AMOEBA method on an elementary reaction step in NanoLuc, an artificial bioluminescence luciferase. We will show how the reaction mechanism is different when calculated in the gas phase, in polarizable continuum medium (PCM), versus in protein AMOEBA models.

Unraveling the Hydration Shell Structure and Dynamics of Group 10 Aqua Ions

We present ab initio simulations based on subsystem DFT of group 10 aqua ions accurately compared against experimental data on hydration structure. Our simulations provide insights into the molecular structures and dynamics of hydration shells, offering recalibrated interpretations of experimental results. We observe a soft, but distinct second hydration shell in Palladium (Pd) due to a balance between thermal fluctuations, metal-water interactions, and hydrogen bonding. Nickel (Ni) and platinum (Pt) exhibit more rigid hydration shells. Notably, our simulations align with experimental findings for Pd, showing axial hydration marked by a broad peak at about 3 Å in the Pd-O radial distribution function, revising the previously sharp "mesoshell" prediction. We introduce the "hydrogen bond dome" concept to describe a resilient network of hydrogen-bonded water molecules around the metal, which plays a critical role in the axial hydration dynamics.

Understanding cellulose pyrolysis via ab initio deep learning potential field

Comprehensive and dynamic studies of cellulose pyrolysis reaction mechanisms are crucial in designing experiments and processes with enhanced safety, efficiency, and sustainability. The details of the pyrolysis mechanism are not readily available from experiments but can be better described via molecular dynamics (MD) simulations. However, the large size of cellulose molecules challenges accurate ab initio MD simulations, while existing reactive force field parameters lack precision. In this work, precise ab initio deep learning potentials field (DPLF) are developed and applied in MD simulations to facilitate the study of cellulose pyrolysis mechanisms. The formation mechanism and production rate of both valuable and greenhouse products from cellulose at temperatures larger than 1073 K are comprehensively described. This study underscores the critical role of advanced simulation techniques, particularly DLPF, in achieving efficient and accurate understanding of cellulose pyrolysis mechanisms, thus promoting wider industrial applications.

Differentiating the Yield of Chemical Reactions Using Parameters in First-Order Kinetic Equations to Identify Elementary Steps That Control the Reactivity from Complicated Reaction Path Networks

The yield of a chemical reaction is obtained by solving its rate equation. This study introduces an approach for differentiating yields by utilizing the parameters of the rate equation, which is expressed as a first-order linear differential equation. The yield derivative for a specific pair of reactants and products is derived by mathematically expressing the rate constant matrix contraction method, which is a simple kinetic analysis method. The parameters of the rate equation are the Gibbs energies of the intermediates and transition states in the reaction path network used to formulate the rate equation. Thus, our approach for differentiating the yield allows a numerical evaluation of the contribution of energy variation to the yield for each intermediate and transition state in the reaction path network. In other words, a comparison of these values automatically extracts the factors affecting the yield from a complicated reaction path network consisting of numerous reaction paths and intermediates. This study verifies the behavior of the proposed approach through numerical experiments on the reaction path networks of a model system and the Rh-catalyzed hydroformylation reaction. Moreover, the possibility of using this approach for designing ligands in organometallic catalysts is discussed.

Exploring Chemical Space Using Ab Initio Hyperreactor Dynamics

In recent years, first-principles exploration of chemical reaction space has provided valuable insights into intricate reaction networks. Here, we introduce ab initio hyperreactor dynamics, which enables rapid screening of the accessible chemical space from a given set of initial molecular species, predicting new synthetic routes that can potentially guide subsequent experimental studies. For this purpose, different hyperdynamics derived bias potentials are applied along with pressure-inducing spherical confinement of the molecular system in ab initio molecular dynamics simulations to efficiently enhance reactivity under mild conditions. To showcase the advantages and flexibility of the hyperreactor approach, we present a systematic study of the method's parameters on a HCN toy model and apply it to a recently introduced experimental model for the prebiotic formation of glycinal and acetamide in interstellar ices, which yields results in line with experimental findings. In addition, we show how the developed framework enables the study of complicated transitions like the first step of a nonenzymatic DNA nucleoside synthesis in an aqueous environment, where the molecular fragmentation problem of earlier nanoreactor approaches is avoided.

Automated Exploration of Reaction Networks and Mechanisms Based on Metadynamics Nanoreactor Simulations

We developed an automated approach to construct a complex reaction network and explore the reaction mechanisms for numerous reactant molecules by integrating several theoretical approaches. Nanoreactor-type molecular dynamics was used to generate possible chemical reactions, in which the metadynamics was used to overcome the reaction barriers, and the semiempirical GFN2-xTB method was used to reduce the computational cost. Reaction events were identified from trajectories using the hidden Markov model based on the evolution of the molecular connectivity. This provided the starting points for further transition-state searches at the electronic structure levels of density functional theory to obtain the reaction mechanism. Finally, the entire reaction network containing multiple pathways was built. The feasibility and efficiency of the automated construction of the reaction network were investigated using the HCHO and NH3 biomolecular reaction and the reaction network for a multispecies system comprising dozens of HCN and H2O molecules. The results indicated that the proposed approach provides a valuable and effective tool for the automated exploration of the reaction networks.

Efficient Acceleration of Reaction Discovery in the Ab Initio Nanoreactor: Phenyl Radical Oxidation Chemistry

Over the years, many computational strategies have been employed to elucidate reaction networks. One of these methods is accelerated molecular dynamics, which can circumvent the expense required in dynamics to find all reactants and products (local minima) and transition states (first-order saddle points) on a potential energy surface (PES) by using fictitious forces that promote reaction events. The ab initio nanoreactor uses these accelerating forces to study large chemical reaction networks from first-principles quantum mechanics. In the initial nanoreactor studies, this acceleration was done through a piston periodic compression potential, which pushes molecules together to induce entropically unfavorable bimolecular reactions. However, the piston is not effective for discovering intramolecular and dissociative reactions, such as those integral to the decomposition channels of phenyl radical oxidation. In fact, the choice of accelerating forces dictates not only the rate of reaction discovery but also the types of reactions discovered; thus, it is critical to understand the biases and efficacies of these forces. In this study, we examine forces using metadynamics, attractive potentials, and local thermostats for accelerating reaction discovery. For each force, we construct a separate phenyl radical combustion reaction network using solely that force in discovery trajectories. We elucidate the enthalpic and entropic trends of each accelerating force and highlight their efficiency in reaction discovery. Comparing the nanoreactor-constructed reaction networks with literature renditions of the phenyl radical combustion PES shows that a combination of accelerating forces is best suited for reaction discovery.

First principles reaction discovery: from the Schrodinger equation to experimental prediction for methane pyrolysis

Our recent success in exploiting graphical processing units (GPUs) to accelerate quantum chemistry computations led to the development of the ab initio nanoreactor, a computational framework for automatic reaction discovery and kinetic model construction. In this work, we apply the ab initio nanoreactor to methane pyrolysis, from automatic reaction discovery to path refinement and kinetic modeling. Elementary reactions occurring during methane pyrolysis are revealed using GPU-accelerated ab initio molecular dynamics simulations. Subsequently, these reaction paths are refined at a higher level of theory with optimized reactant, product, and transition state geometries. Reaction rate coefficients are calculated by transition state theory based on the optimized reaction paths. The discovered reactions lead to a kinetic model with 53 species and 134 reactions, which is validated against experimental data and simulations using literature kinetic models. We highlight the advantage of leveraging local brute force and Monte Carlo sensitivity analysis approaches for efficient identification of important reactions. Both sensitivity approaches can further improve the accuracy of the methane pyrolysis kinetic model. The results in this work demonstrate the power of the ab initio nanoreactor framework for computationally affordable systematic reaction discovery and accurate kinetic modeling.

Further cautionary tales on thermostatting in molecular dynamics: Energy equipartitioning and non-equilibrium processes in gas-phase simulations

Molecular dynamics (MD) simulations of gas-phase chemical reactions are typically carried out on a small number of molecules near thermal equilibrium by means of various thermostatting algorithms. Correct equipartitioning of kinetic energy among translations, rotations, and vibrations of the simulated reactants is critical for many processes occurring in the gas phase. As thermalizing collisions are infrequent in gas-phase simulations, the thermostat has to efficiently reach equipartitioning in the system during equilibration and maintain it throughout the actual simulation. Furthermore, in non-equilibrium simulations where heat is released locally, the action of the thermostat should not lead to unphysical changes in the overall dynamics of the system. Here, we explore issues related to both obtaining and maintaining thermal equilibrium in MD simulations of an exemplary ion-molecule dimerization reaction. We first compare the efficiency of global (Nosé-Hoover and Canonical Sampling through Velocity Rescaling) and local (Langevin) thermostats for equilibrating a system of flexible compounds and find that of these three only the Langevin thermostat achieves equipartition in a reasonable simulation time. We then study the effect of the unphysical removal of latent heat released during simulations involving multiple dimerization events. As the Langevin thermostat does not produce the correct dynamics in the free molecular regime, we only consider the commonly used Nosé-Hoover thermostat, which is shown to effectively cool down the reactants, leading to an overestimation of the dimerization rate. Our findings underscore the importance of thermostatting for the proper thermal initialization of gas-phase systems and the consequences of global thermostatting in non-equilibrium simulations.

Species-selective nanoreactor molecular dynamics simulations based on linear-scaling tight-binding quantum chemical calculations

Here, extensions to quantum chemical nanoreactor molecular dynamics simulations for discovering complex reactive events are presented. The species-selective algorithm, where the nanoreactor effectively works for the selected desired reactants, was introduced to the original scheme. Moreover, for efficient simulations of large model systems with the modified approach, the divide-and-conquer linear-scaling density functional tight-binding method was exploited. Two illustrative applications of the polymerization of propylene and cyclopropane mixtures and the aggregation of sodium chloride from aqueous solutions indicate that species-selective quantum chemical nanoreactor molecular dynamics is a promising method to accelerate the sampling of multicomponent chemical processes proceeding under relatively mild conditions.

Exploring the Quantum Chemical Energy Landscape with GNN-Guided Artificial Force

Artificial force has been proven useful to get over energy barriers and quickly search a large portion of the energy landscape. This work proposes a method based on graph neural networks to optimize the choice of transformation patterns to examine and accelerate energy landscape exploration. In open search from glutathione, the search efficiency was largely improved in comparison to random selection. We also applied transfer learning from glutathione to tuftsin, resulting in further efficiency gains.

A Reactive Molecular Dynamics Study of Chlorinated Organic Compounds. Part I: Force Field Development

This work presents a novel parametrization for the ReaxFF formalism as a means to investigate reaction processes of chlorinated organic compounds. Force field parameters cover the chemical elements C, H, O, Cl and were obtained using a novel optimization approach involving relaxed potential energy surface scans as training targets. The resulting ReaxFF parametrization shows good transferability, as demonstrated on two independent ab initio validation sets. While this first part of our two-paper series focuses on force field parametrization, we apply our parameters to the simulation of chlorinated dibenzofuran formation and decomposition processes in Part II.

Prediction of High-Yielding Single-Step or Cascade Pericyclic Reactions for the Synthesis of Complex Synthetic Targets

Pericyclic reactions, which involve cyclic concerted transition states without ionic or radical intermediates, have been extensively studied since their definition in the 1960s, and the famous Woodward-Hoffmann rules predict their stereoselectivity and chemoselectivity. Here, we describe the application of a fully automated reaction-path search method, that is, the artificial force induced reaction (AFIR), to trace an input compound back to reasonable starting materials through thermally allowed pericyclic reactions via product-based quantum-chemistry-aided retrosynthetic analysis (QCaRA) without using any a priori experimental knowledge. All categories of pericyclic reactions, including cycloadditions, ene reactions, group-transfer, cheletropic, electrocyclic, and sigmatropic reactions, were successfully traced back via concerted reaction pathways, and starting materials were computationally obtained with the correct stereochemistry. Furthermore, AFIR was used to predict whether the identified reaction pathway can be expected to occur in good yield relative to other possible reactions of the identified starting material. In order to showcase its practical utility, this state-of-the-art technology was also applied to the retrosynthetic analysis of a natural product with a relatively high number of atoms (52 atoms: endiandric acid C methyl ester), which was first synthesized by Nicolaou in 1982 and provided the corresponding starting polyenes with the correct stereospecificity via three pericyclic reaction cascades (one Diels-Alder reaction as well as 6π and 8π electrocyclic reactions). Moreover, not only systems that obey the Woodward-Hoffmann rules but also systems that violate these rules, such as those recently calculated by Houk, can be retrosynthesized accurately.

A Reactive Molecular Dynamics Study of Chlorinated Organic Compounds. Part II: A ChemTraYzer Study of Chlorinated Dibenzofuran Formation and Decomposition Processes

In our two-paper series, we first present the development of ReaxFF CHOCl parameters using the recently published ParAMS parametrization tool. In this second part, we update the reactive Molecular Dynamics - Quantum Mechanics coupling scheme ChemTraYzer and combine it with our new ReaxFF parameters from Part I to study formation and decomposition processes of chlorinated dibenzofurans. We introduce a self-learning method for recovering failed transition-state searches that improves the overall ChemTraYzer transition-state search success rate by 10 percentage points to a total of 48 %. With ChemTraYzer, we automatically find and quantify more than 500 reactions using transition state theory and DFT. Among the discovered chlorinated dibenzofuran reactions are numerous reactions that are new to the literature. In three case studies, we discuss the set of reactions that are most relevant to the dibenzofuran literature: (i) bimolecular reactions of the chlorinated-dibenzofuran precursors phenoxy radical and 1,3,5-trichlorobenzene, (ii) dibenzofuran chlorination and pyrolysis, and (iii) oxidation of chlorinated dibenzofurans.

Fully Automated Generation of Prebiotically Relevant Reaction Networks from Optimized Nanoreactor Simulations

The nanoreactor approach first introduced by the group of Martı́nez [Wang et al. Nat. Chem. 2014, 6, 1044-1048] has recently attracted much attention because of its ability to accelerate the discovery of reaction pathways. Here, we provide a comprehensive study of various simulation parameters and present an alternative implementation for the reactivity-enhancing spherical constraint function, as well as for the detection of reaction events. In this context, a fully automated postsimulation evaluation procedure based on RDKit and NetworkX analysis is introduced. The chemical and physical robustness of the procedure is examined by investigating the reactivity of selected homogeneous systems. The optimized procedure is applied at the GFN2-xTB level of theory to a system composed of HCN molecules and argon atoms, acting as a buffer, yielding prebiotically plausible primary and secondary precursors for the synthesis of RNA. Furthermore, the formose reaction network is explored leading to numerous sugar precursors. The discovered compounds reflect experimental findings; however, new synthetic routes and a large collection of exotic, highly reactive molecules are observed, highlighting the predictive power of the nanoreactor approach for unraveling the reactive manifold.

An automated method for graph-based chemical space exploration and transition state finding

Algorithms that automatically explore the chemical space have been limited to chemical systems with a low number of atoms due to expensive involved quantum calculations and the large amount of possible reaction pathways. The method described here presents a novel solution to the problem of chemical exploration by generating reaction networks with heuristics based on chemical theory. First, a second version of the reaction network is determined through molecular graph transformations acting upon functional groups of the reacting. Only transformations that break two chemical bonds and form two new ones are considered, leading to a significant performance enhancement compared to previously presented algorithm. Second, energy barriers for this reaction network are estimated through quantum chemical calculations by a growing string method, which can also identify non-octet species missed during the previous step and further define the reaction network. The proposed algorithm has been successfully applied to five different chemical reactions, in all cases identifying the most important reaction pathways.

Enhanced Sampling Aided Design of Molecular Photoswitches

Advances in the evolving field of atomistic simulations promise important insights for the design and fundamental understanding of novel molecular photoswitches. Here, we use state-of-the-art enhanced simulation techniques to unravel the complex, multistep chemistry of donor-acceptor Stenhouse adducts (DASAs). Our reaction discovery workflow consists of enhanced sampling for efficient chemical space exploration, refinement of newly observed pathways with more accurate ab initio electronic structure calculations, and structural modifications to introduce design principles within future generations of DASAs. We showcase our discovery workflow by not only recovering the full photoswitching mechanism of DASA but also predicting a plethora of new plausible thermal pathways and suggesting a way for their experimental validation. Furthermore, we illustrate the tunability of these newly discovered reactions, leading to a potential avenue for controlling DASA dynamics through multiple external stimuli. Overall, these insights could offer alternative routes to increase the efficiency and control of DASA's photoswitching mechanism, providing new elements to design more complex light-responsive materials.

Self-Supervised Molecular Pretraining Strategy for Low-Resource Reaction Prediction Scenarios

In the face of low-resource reaction training samples, we construct a chemical platform for addressing small-scale reaction prediction problems. Using a self-supervised pretraining strategy called MAsked Sequence to Sequence (MASS), the Transformer model can absorb the chemical information of about 1 billion molecules and then fine-tune on a small-scale reaction prediction. To further strengthen the predictive performance of our model, we combine MASS with the reaction transfer learning strategy. Here, we show that the average improved accuracies of the Transformer model can reach 14.07, 24.26, 40.31, and 57.69% in predicting the Baeyer-Villiger, Heck, C-C bond formation, and functional group interconversion reaction data sets, respectively, marking an important step to low-resource reaction prediction.

O-Acetyl Migration within the Sialic Acid Side Chain: A Mechanistic Study Using the Ab Initio Nanoreactor

Many disease-causing viruses target sialic acids on the surface of host cells. Some viruses bind preferentially to sialic acids with O-acetyl modification at the hydroxyl group of C7, C8, or C9 on the glycerol-like side chain. Studies of proteins binding to sialosides containing O-acetylated sialic acids are crucial in understanding the related diseases but experimentally difficult due to the lability of the ester group. We recently showed that O-acetyl migration among hydroxyl groups of C7, C8, and C9 in sialic acids occurs in all directions in a pH-dependent manner. In the current study, we elucidate a full mechanistic pathway for the migration of O-acetyl among C7, C8, and C9. We used an ab initio nanoreactor to explore potential reaction pathways and density functional theory, pK_a calculations, and umbrella sampling to investigate elementary steps of interest. We found that when a base is present, migration is easy in any direction and involves three key steps: deprotonation of the hydroxyl group, cyclization between the two carbons, and the migration of the O-acetyl group. This dynamic equilibrium may play a defensive role against pathogens that evolve to gain entry to the cell by binding selectively to one acetylation state.

State-averaged CASSCF with polarizable continuum model for studying photoreactions in solvents: Energies, analytical nuclear gradients, and non-adiabatic couplings

This paper presents state-averaged complete active space self-consistent field in polarizable continuum model (PCM) for studies of photoreactions in solvents. The wavefunctions of the solute and the PCM surface charges of the solvent are optimized simultaneously such that the state-averaged free energy is variationally minimized. The method supports both fixed weights and dynamic weights where the weights are automatically adjusted based on the energy gaps. The corresponding analytical nuclear gradients and non-adiabatic couplings are also derived. Furthermore, we show how the new method can be entirely formulated in terms of seven basic operations, which allows the implementation to benefit from existing high-performance libraries on graphical processing units. Results demonstrating the accuracy and performance of the implementation are presented and discussed. We also apply the new method to the study of minimal conical intersection search and photoreaction energy pathways in solvents. Effects from the polarity of the solvents and different formulas of dynamic weights are compared and discussed.

Discover, Sample, and Refine: Exploring Chemistry with Enhanced Sampling Techniques

Over the last few decades, enhanced sampling methods have been continuously improved. Here, we exploit this progress and propose a modular workflow for blind reaction discovery and determination of reaction paths. In a three-step strategy, at first we use a collective variable derived from spectral graph theory in conjunction with the explore variant of the on-the-fly probability enhanced sampling method to drive reaction discovery runs. Once different chemical products are determined, we construct an ad-hoc neural network-based collective variable to improve sampling, and finally we refine the results using the free energy perturbation theory and a more accurate Hamiltonian. We apply this strategy to both intramolecular and intermolecular reactions. Our workflow requires minimal user input and extends the power of ab initio molecular dynamics to explore and characterize the reaction space.

Efficient Reaction Space Exploration with ChemTraYzer-TAD

The development of a reaction model is often a time-consuming process, especially if unknown reactions have to be found and quantified. To alleviate the reaction modeling process, automated procedures for reaction space exploration are highly desired. We present ChemTraYzer-TAD, a new reactive molecular dynamics acceleration technique aimed at efficient reaction space exploration. The new method is based on the basin confinement strategy known from the temperature-accelerated dynamics (TAD) acceleration method. Our method features integrated ChemTraYzer bond-order processing steps for the automatic and on-the-fly determination of the positions of virtual walls in configuration space that confine the system in a potential energy basin. We use the example of 1,3-dioxolane-4-hydroperoxide-2-yl radical oxidation to show that ChemTraYzer-TAD finds more than 100 different parallel reactions for the given set of reactants in less than 2 ns of simulation time. Among the many observed reactions, ChemTraYzer-TAD finds the expected typical low-temperature reactions despite the use of extremely high simulation temperatures up to 5000 K. Our method also finds a new concerted β-scission plus O₂ addition with a lower reaction barrier than the literature-known and so-far dominant β-scission.

CARNOT: a Fragment-Based Direct Molecular Dynamics and Virtual-Reality Simulation Package for Reactive Systems

Traditionally, the study of reaction mechanisms of complex reaction systems such as combustion has been performed on an individual basis by optimizations of transition structure and minimum energy path or by reaction dynamics trajectory calculations for one elementary reaction at a time. It is effective, but time-consuming, whereas important and unexpected processes could have been missed. In this article, we present a direct molecular dynamics (DMD) approach and a virtual-reality simulation program, CARNOT, in which plausible chemical reactions are simulated simultaneously at finite temperature and pressure conditions. A key concept of the present ab initio molecular dynamics method is to partition a large, chemically reactive system into molecular fragments that can be adjusted on the fly of a DMD simulation. The theory represents an extension of the explicit polarization method to reactive events, called ReX-Pol. We propose a highest-and-lowest adapted-spin approximation to define the local spins of individual fragments, rather than treating the entire system by a delocalized wave function. Consequently, the present ab initio DMD can be applied to reactive systems consisting of an arbitrarily varying number of closed and open-shell fragments such as free radicals, zwitterions, and separate ions found in combustion and other reactions. A graph-data structure algorithm was incorporated in CARNOT for the analysis of reaction networks, suitable for reaction mechanism reduction. Employing the PW91 density functional theory and the 6-31+G(d) basis set, the capabilities of the CARNOT program were illustrated by a combustion reaction, consisting of 28 650 atoms, and by reaction network analysis that revealed a range of mechanistic and dynamical events. The method may be useful for applications to other types of complex reactions.

New Tools for Taming Complex Reaction Networks: The Unimolecular Decomposition of Indole Revisited

The level of detail attained in the computational description of reaction mechanisms can be vastly improved through tools for automated chemical space exploration, particularly for systems of small to medium size. Under this approach, the unimolecular decomposition landscape for indole was explored through the automated reaction mechanism discovery program AutoMeKin. Nevertheless, the sheer complexity of the obtained mechanisms might be a hindrance regarding their chemical interpretation. In this spirit, the new Python library amk-tools has been designed to read and manipulate complex reaction networks, greatly simplifying their overall analysis. The package provides interactive dashboards featuring visualizations of the network, the three-dimensional (3D) molecular structures and vibrational normal modes of all chemical species, and the corresponding energy profiles for selected pathways. The combination of the joined mechanism generation and postprocessing workflow with the rich chemistry of indole decomposition enabled us to find new details of the reaction (obtained at the CCSD(T)/aug-cc-pVTZ//M06-2X/MG3S level of theory) that were not reported before: (i) 16 pathways leading to the formation of HCN and NH₃ (via amino radical); (ii) a barrierless reaction between methylene radical and phenyl isocyanide, which might be an operative mechanism under the conditions of the interstellar medium; and (iii) reaction channels leading to both hydrogen cyanide and hydrogen isocyanide, of potential astrochemical interest as the computed HNC/HCN ratios greatly exceed the calculated equilibrium value at very low temperatures. The reported reaction networks can be very valuable to supplement databases of kinetic data, which is of remarkable interest for pyrolysis and astrochemical studies.

Computational Discovery of Transition-metal Complexes: From High-throughput Screening to Machine Learning

Transition-metal complexes are attractive targets for the design of catalysts and functional materials. The behavior of the metal-organic bond, while very tunable for achieving target properties, is challenging to predict and necessitates searching a wide and complex space to identify needles in haystacks for target applications. This review will focus on the techniques that make high-throughput search of transition-metal chemical space feasible for the discovery of complexes with desirable properties. The review will cover the development, promise, and limitations of "traditional" computational chemistry (i.e., force field, semiempirical, and density functional theory methods) as it pertains to data generation for inorganic molecular discovery. The review will also discuss the opportunities and limitations in leveraging experimental data sources. We will focus on how advances in statistical modeling, artificial intelligence, multiobjective optimization, and automation accelerate discovery of lead compounds and design rules. The overall objective of this review is to showcase how bringing together advances from diverse areas of computational chemistry and computer science have enabled the rapid uncovering of structure-property relationships in transition-metal chemistry. We aim to highlight how unique considerations in motifs of metal-organic bonding (e.g., variable spin and oxidation state, and bonding strength/nature) set them and their discovery apart from more commonly considered organic molecules. We will also highlight how uncertainty and relative data scarcity in transition-metal chemistry motivate specific developments in machine learning representations, model training, and in computational chemistry. Finally, we will conclude with an outlook of areas of opportunity for the accelerated discovery of transition-metal complexes.

AutoMeKin2021: An open-source program for automated reaction discovery

AutoMeKin2021 is an updated version of tsscds2018, a program for the automated discovery of reaction mechanisms (J. Comput. Chem. 2018, 39, 1922). This release features a number of new capabilities: rare-event molecular dynamics simulations to enhance reaction discovery, extension of the original search algorithm to study van der Waals complexes, use of chemical knowledge, a new search algorithm based on bond-order time series analysis, statistics of the chemical reaction networks, a web application to submit jobs, and other features. The source code, manual, installation instructions and the website link are available at: https://rxnkin.usc.es/index.php/AutoMeKin.

Machine Learning for Chemical Reactions

Machine learning (ML) techniques applied to chemical reactions have a long history. The present contribution discusses applications ranging from small molecule reaction dynamics to computational platforms for reaction planning. ML-based techniques can be particularly relevant for problems involving both computation and experiments. For one, Bayesian inference is a powerful approach to develop models consistent with knowledge from experiments. Second, ML-based methods can also be used to handle problems that are formally intractable using conventional approaches, such as exhaustive characterization of state-to-state information in reactive collisions. Finally, the explicit simulation of reactive networks as they occur in combustion has become possible using machine-learned neural network potentials. This review provides an overview of the questions that can and have been addressed using machine learning techniques, and an outlook discusses challenges in this diverse and stimulating field. It is concluded that ML applied to chemistry problems as practiced and conceived today has the potential to transform the way with which the field approaches problems involving chemical reactions, in both research and academic teaching.

Putting Density Functional Theory to the Test in Machine-Learning-Accelerated Materials Discovery

Accelerated discovery with machine learning (ML) has begun to provide the advances in efficiency needed to overcome the combinatorial challenge of computational materials design. Nevertheless, ML-accelerated discovery both inherits the biases of training data derived from density functional theory (DFT) and leads to many attempted calculations that are doomed to fail. Many compelling functional materials and catalytic processes involve strained chemical bonds, open-shell radicals and diradicals, or metal-organic bonds to open-shell transition-metal centers. Although promising targets, these materials present unique challenges for electronic structure methods and combinatorial challenges for their discovery. In this Perspective, we describe the advances needed in accuracy, efficiency, and approach beyond what is typical in conventional DFT-based ML workflows. These challenges have begun to be addressed through ML models trained to predict the results of multiple methods or the differences between them, enabling quantitative sensitivity analysis. For DFT to be trusted for a given data point in a high-throughput screen, it must pass a series of tests. ML models that predict the likelihood of calculation success and detect the presence of strong correlation will enable rapid diagnoses and adaptation strategies. These "decision engines" represent the first steps toward autonomous workflows that avoid the need for expert determination of the robustness of DFT-based materials discoveries.

The non-adiabatic nanoreactor: towards the automated discovery of photochemistry

The ab initio nanoreactor has previously been introduced to automate reaction discovery for ground state chemistry. In this work, we present the nonadiabatic nanoreactor, an analogous framework for excited state reaction discovery. We automate the study of nonadiabatic decay mechanisms of molecules by probing the intersection seam between adiabatic electronic states with hyper-real metadynamics, sampling the branching plane for relevant conical intersections, and performing seam-constrained path searches. We illustrate the effectiveness of the nonadiabatic nanoreactor by applying it to benzene, a molecule with rich photochemistry and a wide array of photochemical products. Our study confirms the existence of several types of S₀/S₁ and S₁/S₂ conical intersections which mediate access to a variety of ground state stationary points. We elucidate the connections between conical intersection energy/topography and the resulting photoproduct distribution, which changes smoothly along seam space segments. The exploration is performed with minimal user input, and the protocol requires no previous knowledge of the photochemical behavior of a target molecule. We demonstrate that the nonadiabatic nanoreactor is a valuable tool for the automated exploration of photochemical reactions and their mechanisms.

Nitromethane Decomposition via Automated Reaction Discovery and an Ab Initio Corrected Kinetic Model

We explore the systematic construction of kinetic models from in silico reaction data for the decomposition of nitromethane. Our models are constructed in a computationally affordable manner by using reactions discovered through accelerated molecular dynamics simulations using the ReaxFF reactive force field. The reaction paths are then optimized to determine reaction rate parameters. We introduce a reaction barrier correction scheme that combines accurate thermochemical data from density functional theory with ReaxFF minimal energy paths. We validate our models across different thermodynamic regimes, showing predictions of gas phase CO and NO concentrations and high-pressure induction times that are similar to experimental data. The kinetic models are analyzed to find fundamental decomposition reactions in different thermodynamic regimes.

Simple position and orientation preconditioning scheme for minimum energy path calculations

Minimum-energy path (MEP) calculations, such as those typified by the nudged elastic band method, require input of reactant and product molecular configurations at initialization. In the case of reactions involving more than one molecule, generating initial reactant and product configurations requires careful consideration of the relative position and orientations of the reactive molecules in order to ensure that the resulting MEP calculation proceeds without converging on an alternative reaction-path, and without requiring excessive numbers of optimization iterations; as such, this initial system set-up is most commonly performed "by hand," with an expert user arranging reactive molecules in space to ensure that the following MEP calculation runs smoothly. In this Article, we introduce a simple preconditioning scheme which replaces this labor-intensive, human-knowledge-based step with an automated deterministic computational scheme. In our approach, initial reactant and product configurations are generated such that steric hindrance between reactive molecules is minimized in the reactant and product configurations, while also simultaneously requiring minimal structural differences between the reactants and products. The method is demonstrated using a benchmark test-set of >3400 organic molecular reactions, where comparison of the reactant/product configurations generated using our approach compare very well to initial configurations which were generated on an ad hoc basis.

Data augmentation and transfer learning strategies for reaction prediction in low chemical data regimes

An effective and rapid deep learning method to predict chemical reactions contributes to the research and development of organic chemistry and drug discovery.

Reduced scaling extended multi-state CASPT2 (XMS-CASPT2) using supporting subspaces and tensor hyper-contraction

We present a reduced scaling formulation of the extended multi-state CASPT2 (XMS-CASPT2) method, which is based on our recently developed state-specific CASPT2 (SS-CASPT2) formulation using supporting subspaces and tensor hyper-contraction. By using these two techniques, the off-diagonal elements of the effective Hamiltonian can be computed with only O(N³) operations and O(N²) memory, where N is the number of basis functions. This limits the overall computational scaling to O(N⁴) operations and O(N²) memory. Thus, excited states can now be obtained at the same reduced (relative to previous algorithms) scaling we achieved for SS-CASPT2. In addition, we also investigate how the energy denominators can be factorized with the Laplace quadrature when some of the denominators are negative, which is critical for excited state calculations. An efficient implementation of the method has been developed using graphical processing units while also exploiting spatial sparsity in tensor operations. We benchmark the accuracy of the new method by comparison to non-THC formulated XMS-CASPT2 for the excited states of various molecules. In our tests, the THC approximation introduces negligible errors (≈0.01 eV) compared to the non-THC reference method. Scaling behavior and computational timings are presented to demonstrate performance. The new method is also interfaced with quantum mechanics/molecular mechanics (QM/MM). In an example study of green fluorescent protein, we show how the XMS-CASPT2 potential energy surfaces and excitation energies are affected by increasing the size of the QM region up to 278 QM atoms with more than 2300 basis functions.