Removing External Degrees of Freedom from Transition-State Search Methods using Quaternions

Highly Accurate Prediction of NMR Chemical Shifts from Low-Level Quantum Mechanics Calculations Using Machine Learning

Theoretical predictions of NMR chemical shifts from first-principles can greatly facilitate experimental interpretation and structure identification of molecules in gas, solution, and solid-state phases. However, accurate prediction of chemical shifts using the gold-standard coupled cluster with singles, doubles, and perturbative triple excitations [CCSD(T)] method with a complete basis set (CBS) can be prohibitively expensive. By contrast, machine learning (ML) methods offer inexpensive alternatives for chemical shift predictions but are hampered by generalization to molecules outside the original training set. Here, we propose several new ideas in ML of the chemical shift prediction for H, C, N, and O that first introduce a novel feature representation, based on the atomic chemical shielding tensors within a molecular environment using an inexpensive quantum mechanics (QM) method, and train it to predict NMR chemical shieldings of a high-level composite theory that approaches the accuracy of CCSD(T)/CBS. In addition, we train the ML model through a new progressive active learning workflow that reduces the total number of expensive high-level composite calculations required while allowing the model to continuously improve on unseen data. Furthermore, the algorithm provides an error estimation, signaling potential unreliability in predictions if the error is large. Finally, we introduce a novel approach to keep the rotational invariance of the features using tensor environment vectors (TEVs) that yields a ML model with the highest accuracy compared to a similar model using data augmentation. We illustrate the predictive capacity of the resulting inexpensive shift machine learning (iShiftML) models across several benchmarks, including unseen molecules in the NS372 data set, gas-phase experimental chemical shifts for small organic molecules, and much larger and more complex natural products in which we can accurately differentiate between subtle diastereomers based on chemical shift assignments.

Improved Initialization of Optimal Path Calculations Using Sequential Traversal over the Image-Dependent Pair Potential Surface

In reaction path optimization, such as the calculation of a minimum energy path (MEP) between given reactant and product configurations of atoms, it is advantageous to start with an initial guess where the close proximity of atoms is avoided and bonds are not unnecessarily broken only to be reformed later. When the configurations of the atoms are described with Cartesian coordinates, a linear interpolation between the end points can be problematic, and a better option is provided by the so-called image dependent pair potential (IDPP) approach where interpolated pairwise distances are generated to form an objective function that can be used to construct an improved initial path. When started with a linear interpolation, this method can, however, still lead to unnecessary bond breaking in, for example, reactions in which a molecular subgroup undergoes significant rotation. In the method presented here, this problem is addressed by constructing the path gradually, introducing images sequentially starting from the vicinity of the end points while the distance between images in the central region is larger. The distribution of images is controlled by systematically scaling the tightness of springs acting between the images until the desired number of images is obtained, and they are evenly spaced. This procedure generates an initial path on the IDPP surface, a task that requires negligible computational effort, as no evaluation of the energy of the system is needed. The calculation of the MEP, typically using electronic structure calculations, is then subsequently carried out in a way that makes efficient use of parallel computing with the nudged elastic band method. Several examples of reactions are given where the linear interpolation IDPP (LI-IDPP) method yields problematic paths with unnecessary bond breaking in some of the intermediate images, while the sequential IDPP (S-IDPP) method yields paths that are significantly closer to realistic MEPs.

Generalized Energy-Based Fragmentation Approach for the Electronic Emission Spectra of Large Systems

The excited-state (ES) geometry optimization and electronic emission (fluorescence and phosphorescence) spectra and the ES vibrational spectra of large systems are great challenges in quantum chemistry. In this work, we develop a generalized energy-based fragmentation (GEBF) approach to compute the localized ES structures and vibrational frequencies of large systems. In this approach, the ES energy derivatives (gradients or Hessians) for a localized ES of a large system can be obtained by combining the ES energy derivatives of the corresponding active subsystems (including local excitation center) and the ground-state energy derivatives of inactive subsystems. Two strategies are adopted to overcome two difficulties from state-classification and state-tracking for treating specific ESs. First, for state-classification, we develop an improved density-based spatial clustering applied with noise algorithm with a modified transition orbital projection (TOP) algorithm, which allow a certain ES energy and energy derivatives of the whole system to be calculated with different ES energies and energy derivatives of active subsystems. Furthermore, we also employ the TOP algorithm for tracking the ESs in their geometry optimizations at the time-dependent density functional theory (TDDFT) level. Then, the GEBF approach is applied to investigate the optimized ES geometries or ES vibrational frequencies for two typical systems. Our results show that the cost-effective GEBF approach can accurately reproduce the TDDFT fluorescence spectra of the cytosine derivative and the experimental phosphorescence spectra of the β-cyclodextrin derivative. The GEBF approach is expected to be routinely applied to investigate the electronic emission spectra of very large systems with local chromophores.

Nudged Elastic Band Method for Molecular Reactions Using Energy-Weighted Springs Combined with Eigenvector Following

The climbing image nudged elastic band method (CI-NEB) is used to identify reaction coordinates and to find saddle points representing transition states of reactions. It can make efficient use of parallel computing as the calculations of the discretization points, the so-called images, can be carried out simultaneously. In typical implementations, the images are distributed evenly along the path by connecting adjacent images with equally stiff springs. However, for systems with a high degree of flexibility, this can lead to poor resolution near the saddle point. By making the spring constants increase with energy, the resolution near the saddle point is improved. To assess the performance of this energy-weighted CI-NEB method, calculations are carried out for a benchmark set of 121 molecular reactions. The performance of the method is analyzed with respect to the input parameters. Energy-weighted springs are found to greatly improve performance and result in successful location of the saddle points in less than a thousand energy and force evaluations on average (about a hundred per image) using the same set of parameter values for all of the reactions. Even better performance is obtained by stopping the calculation before full convergence and complete the saddle point search using an eigenvector following method starting from the location of the climbing image. This combination of methods, referred to as NEB-TS, turns out to be robust and highly efficient as it reduces the average number of energy and force evaluations down to a third, to 305. An efficient and flexible implementation of these methods has been made available in the ORCA software.

Conformational analysis of tannic acid: Environment effects in electronic and reactivity properties

Polyphenols are natural molecules of crucial importance in many applications, of which tannic acid (TA) is one of the most abundant and established. Most high-value applications require precise control of TA interactions with the system of interest. However, the molecular structure of TA is still not comprehended at the atomic level, of which all electronic and reactivity properties depend. Here, we combine an enhanced sampling global optimization method with density functional theory (DFT)-based calculations to explore the conformational space of TA assisted by unsupervised machine learning visualization and then investigate its lowest energy conformers. We study the external environment's effect on the TA structure and properties. We find that vacuum favors compact structures by stabilizing peripheral atoms' weak interactions, while in water, the molecule adopts more open conformations. The frontier molecular orbitals of the conformers with the lowest harmonic vibrational free energy have a HOMO-LUMO energy gap of 2.21 (3.27) eV, increasing to 2.82 (3.88) eV in water, at the DFT generalized gradient approximation (and hybrid) level of theory. Structural differences also change the distribution of potential reactive sites. We establish the fundamental importance of accurate structural consideration in determining TA and related polyphenol interactions in relevant technological applications.

Simultaneously improving reaction coverage and computational cost in automated reaction prediction tasks

Automated reaction prediction has the potential to elucidate complex reaction networks for applications ranging from combustion to materials degradation, but computational cost and inconsistent reaction coverage are still obstacles to exploring deep reaction networks. Here we show that cost can be reduced and reaction coverage can be increased simultaneously by relatively straightforward modifications of the reaction enumeration, geometry initialization and transition state convergence algorithms that are common to many prediction methodologies. These components are implemented in the context of yet another reaction program (YARP), our reaction prediction package with which we report reaction discovery benchmarks for organic single-step reactions, thermal degradation of a γ-ketohydroperoxide, and competing ring-closures in a large organic molecule. Compared with recent benchmarks, YARP (re)discovers both established and unreported reaction pathways and products while simultaneously reducing the cost of reaction characterization by nearly 100-fold and increasing convergence of transition states. This combination of ultra-low cost and high reaction coverage creates opportunities to explore the reactivity of larger systems and more complex reaction networks for applications such as chemical degradation, where computational cost is a bottleneck.

QM7-X, a comprehensive dataset of quantum-mechanical properties spanning the chemical space of small organic molecules

We introduce QM7-X, a comprehensive dataset of 42 physicochemical properties for ≈4.2 million equilibrium and non-equilibrium structures of small organic molecules with up to seven non-hydrogen (C, N, O, S, Cl) atoms. To span this fundamentally important region of chemical compound space (CCS), QM7-X includes an exhaustive sampling of (meta-)stable equilibrium structures-comprised of constitutional/structural isomers and stereoisomers, e.g., enantiomers and diastereomers (including cis-/trans- and conformational isomers)-as well as 100 non-equilibrium structural variations thereof to reach a total of ≈4.2 million molecular structures. Computed at the tightly converged quantum-mechanical PBE0+MBD level of theory, QM7-X contains global (molecular) and local (atom-in-a-molecule) properties ranging from ground state quantities (such as atomization energies and dipole moments) to response quantities (such as polarizability tensors and dispersion coefficients). By providing a systematic, extensive, and tightly-converged dataset of quantum-mechanically computed physicochemical properties, we expect that QM7-X will play a critical role in the development of next-generation machine-learning based models for exploring greater swaths of CCS and performing in silico design of molecules with targeted properties.

Globally Optimal Catalytic Fields - Inverse Design of Abstract Embeddings for Maximum Reaction Rate Acceleration

The search for, and understanding of, good catalysts for chemical reactions is a central issue for chemists. Here, we present first steps toward developing a general computational framework to better support this task. This framework combines efficient, unbiased global optimization techniques with an abstract representation of the catalytic environment, to shrink the search space. To analyze the resulting catalytic embeddings, we employ dimensionality reduction and clustering techniques. This not only provides an inverse design approach to new catalytic embeddings but also illuminates the actual interactions behind catalytic effects. All this is illustrated here with a strictly electrostatic model for the environment and with two versions of a selected example reaction. We close with detailed discussions of future improvements of our framework.

Automated Transition State Search and Its Application to Diverse Types of Organic Reactions

Transition state search is at the center of multiple types of computational chemical predictions related to mechanistic investigations, reactivity and regioselectivity predictions, and catalyst design. The process of finding transition states in practice is, however, a laborious multistep operation that requires significant user involvement. Here, we report a highly automated workflow designed to locate transition states for a given elementary reaction with minimal setup overhead. The only essential inputs required from the user are the structures of the separated reactants and products. The seamless workflow combining computational technologies from the fields of cheminformatics, molecular mechanics, and quantum chemistry automatically finds the most probable correspondence between the atoms in the reactants and the products, generates a transition state guess, launches a transition state search through a combined approach involving the relaxing string method and the quadratic synchronous transit, and finally validates the transition state via the analysis of the reactive chemical bonds and imaginary vibrational frequencies as well as by the intrinsic reaction coordinate method. Our approach does not target any specific reaction type, nor does it depend on training data; instead, it is meant to be of general applicability for a wide variety of reaction types. The workflow is highly flexible, permitting modifications such as a choice of accuracy, level of theory, basis set, or solvation treatment. Successfully located transition states can be used for setting up transition state guesses in related reactions, saving computational time and increasing the probability of success. The utility and performance of the method are demonstrated in applications to transition state searches in reactions typical for organic chemistry, medicinal chemistry, and homogeneous catalysis research. In particular, applications of our code to Michael additions, hydrogen abstractions, Diels-Alder cycloadditions, carbene insertions, and an enzyme reaction model involving a molybdenum complex are shown and discussed.

The atomic simulation environment-a Python library for working with atoms

The atomic simulation environment (ASE) is a software package written in the Python programming language with the aim of setting up, steering, and analyzing atomistic simulations. In ASE, tasks are fully scripted in Python. The powerful syntax of Python combined with the NumPy array library make it possible to perform very complex simulation tasks. For example, a sequence of calculations may be performed with the use of a simple 'for-loop' construction. Calculations of energy, forces, stresses and other quantities are performed through interfaces to many external electronic structure codes or force fields using a uniform interface. On top of this calculator interface, ASE provides modules for performing many standard simulation tasks such as structure optimization, molecular dynamics, handling of constraints and performing nudged elastic band calculations.

Path optimization by a variational reaction coordinate method. I. Development of formalism and algorithms

The development of algorithms to optimize reaction pathways between reactants and products is an active area of study. Existing algorithms typically describe the path as a discrete series of images (chain of states) which are moved downhill toward the path, using various reparameterization schemes, constraints, or fictitious forces to maintain a uniform description of the reaction path. The Variational Reaction Coordinate (VRC) method is a novel approach that finds the reaction path by minimizing the variational reaction energy (VRE) of Quapp and Bofill. The VRE is the line integral of the gradient norm along a path between reactants and products and minimization of VRE has been shown to yield the steepest descent reaction path. In the VRC method, we represent the reaction path by a linear expansion in a set of continuous basis functions and find the optimized path by minimizing the VRE with respect to the linear expansion coefficients. Improved convergence is obtained by applying constraints to the spacing of the basis functions and coupling the minimization of the VRE to the minimization of one or more points along the path that correspond to intermediates and transition states. The VRC method is demonstrated by optimizing the reaction path for the Müller-Brown surface and by finding a reaction path passing through 5 transition states and 4 intermediates for a 10 atom Lennard-Jones cluster.