Predicting Long-Time-Scale Kinetics under Variable Experimental Conditions with Kinetica.jl.

Predicting the degradation processes of molecules over long time scales is a key aspect of industrial materials design. However, it is made computationally challenging by the need to construct large networks of chemical reactions that are relevant to the experimental conditions that kinetic models must mirror, with every reaction requiring accurate kinetic data. Here, we showcase Kinetica.jl, a new software package for constructing large-scale chemical reaction networks in a fully automated fashion by exploring chemical reaction space with a kinetics-driven algorithm; coupled to efficient machine-learning models of activation energies for sampled elementary reactions, we show how this approach readily enables generation and kinetic characterization of networks containing ∼103 chemical species and ≃104-105 reactions. Symbolic-numeric modeling of the generated reaction networks is used to allow for flexible, efficient computation of kinetic profiles under experimentally realizable conditions such as continuously variable temperature regimes, enabling direct connection between bottom-up reaction networks and experimental observations. Highly efficient propagation of long-time-scale kinetic profiles is required for automated reaction network refinement and is enabled here by a new discrete kinetic approximation. The resulting Kinetica.jl simulation package therefore enables automated generation, characterization, and long-time-scale modeling of complex chemical reaction systems. We demonstrate this for hydrocarbon pyrolysis simulated over time scales of seconds, using transient temperature profiles representing those of tubular flow reactor experiments.

Machine Learning Interatomic Potentials for Reactive Hydrogen Dynamics at Metal Surfaces Based on Iterative Refinement of Reaction Probabilities.

The reactive chemistry of molecular hydrogen at surfaces, notably dissociative sticking and hydrogen evolution, plays a crucial role in energy storage and fuel cells. Theoretical studies can help to decipher underlying mechanisms and reaction design, but studying dynamics at surfaces is computationally challenging due to the complex electronic structure at interfaces and the high sensitivity of dynamics to reaction barriers. In addition, ab initio molecular dynamics, based on density functional theory, is too computationally demanding to accurately predict reactive sticking or desorption probabilities, as it requires averaging over tens of thousands of initial conditions. High-dimensional machine learning-based interatomic potentials are starting to be more commonly used in gas-surface dynamics, yet robust approaches to generate reliable training data and assess how model uncertainty affects the prediction of dynamic observables are not well established. Here, we employ ensemble learning to adaptively generate training data while assessing model performance with full uncertainty quantification (UQ) for reaction probabilities of hydrogen scattering on different copper facets. We use this approach to investigate the performance of two message-passing neural networks, SchNet and PaiNN. Ensemble-based UQ and iterative refinement allow us to expose the shortcomings of the invariant pairwise-distance-based feature representation in the SchNet model for gas-surface dynamics.

Assessing Mixed Quantum-Classical Molecular Dynamics Methods for Nonadiabatic Dynamics of Molecules on Metal Surfaces.

Mixed quantum-classical (MQC) methods for simulating the dynamics of molecules at metal surfaces have the potential to accurately and efficiently provide mechanistic insight into reactive processes. Here, we introduce simple two-dimensional models for the scattering of diatomic molecules at metal surfaces based on recently published electronic structure data. We apply several MQC methods to investigate their ability to capture how nonadiabatic effects influence molecule-metal energy transfer during the scattering process. Specifically, we compare molecular dynamics with electronic friction, Ehrenfest dynamics, independent electron surface hopping, and the broadened classical master equation approach. In the case of independent electron surface hopping, we implement a simple decoherence correction approach and assess its impact on vibrationally inelastic scattering. Our results show that simple, low-dimensional models can be used to qualitatively capture experimentally observed vibrational energy transfer and provide insight into the relative performance of different MQC schemes. We observe that all approaches predict similar kinetic energy dependence but return different vibrational energy distributions. Finally, by varying the molecule-metal coupling, we can assess the coupling regime in which some MQC methods become unsuitable.

Generating Protein Folding Trajectories Using Contact-Map-Driven Directed Walks.

Recent advances in machine learning methods have had a significant impact on protein structure prediction, but accurate generation and characterization of protein-folding pathways remains intractable. Here, we demonstrate how protein folding trajectories can be generated using a directed walk strategy operating in the space defined by the residue-level contact-map. This double-ended strategy views protein folding as a series of discrete transitions between connected minima on the potential energy surface. Subsequent reaction-path analysis for each transition enables thermodynamic and kinetic characterization of each protein-folding path. We validate the protein-folding paths generated by our discretized-walk strategy against direct molecular dynamics simulations for a series of model coarse-grained proteins constructed from hydrophobic and polar residues. This comparison demonstrates that ranking discretized paths based on the intermediate energy barriers provides a convenient route to identifying physically sensible folding ensembles. Importantly, by using directed walks in the protein contact-map space, we circumvent several of the traditional challenges associated with protein-folding studies, namely, long time scales required and the choice of a specific order parameter to drive the folding process. As such, our approach offers a useful new route for studying the protein-folding problem.

Efficient implementation and performance analysis of the independent electron surface hopping method for dynamics at metal surfaces.

Independent electron surface hopping (IESH) is a computational algorithm for simulating the mixed quantum-classical molecular dynamics of adsorbate atoms and molecules interacting with metal surfaces. It is capable of modeling the nonadiabatic effects of electron-hole pair excitations on molecular dynamics. Here, we present a transparent, reliable, and efficient implementation of IESH, demonstrating its ability to predict scattering and desorption probabilities across a variety of systems, ranging from model Hamiltonians to full dimensional atomistic systems. We further show how the algorithm can be modified to account for the application of an external bias potential, comparing its accuracy to results obtained using the hierarchical quantum master equation. Our results show that IESH is a practical method for modeling coupled electron-nuclear dynamics at metal surfaces, especially for highly energetic scattering events.

Graph-Driven Reaction Discovery: Progress, Challenges, and Future Opportunities.

Graph-based descriptors, such as bond-order matrices and adjacency matrices, offer a simple and compact way of categorizing molecular structures; furthermore, such descriptors can be readily used to catalog chemical reactions (i.e., bond-making and -breaking). As such, a number of graph-based methodologies have been developed with the goal of automating the process of generating chemical reaction network models describing the possible mechanistic chemistry in a given set of reactant species. Here, we outline the evolution of these graph-based reaction discovery schemes, with particular emphasis on more recent methods incorporating graph-based methods with semiempirical and ab initio electronic structure calculations, minimum-energy path refinements, and transition state searches. Using representative examples from homogeneous catalysis and interstellar chemistry, we highlight how these schemes increasingly act as "virtual reaction vessels" for interrogating mechanistic questions. Finally, we highlight where challenges remain, including issues of chemical accuracy and calculation speeds, as well as the inherent challenge of dealing with the vast size of accessible chemical reaction space.

NQCDynamics.jl: A Julia package for nonadiabatic quantum classical molecular dynamics in the condensed phase.

Accurate and efficient methods to simulate nonadiabatic and quantum nuclear effects in high-dimensional and dissipative systems are crucial for the prediction of chemical dynamics in the condensed phase. To facilitate effective development, code sharing, and uptake of newly developed dynamics methods, it is important that software implementations can be easily accessed and built upon. Using the Julia programming language, we have developed the NQCDynamics.jl package, which provides a framework for established and emerging methods for performing semiclassical and mixed quantum-classical dynamics in the condensed phase. The code provides several interfaces to existing atomistic simulation frameworks, electronic structure codes, and machine learning representations. In addition to the existing methods, the package provides infrastructure for developing and deploying new dynamics methods, which we hope will benefit reproducibility and code sharing in the field of condensed phase quantum dynamics. Herein, we present our code design choices and the specific Julia programming features from which they benefit. We further demonstrate the capabilities of the package on two examples of chemical dynamics in the condensed phase: the population dynamics of the spin-boson model as described by a wide variety of semiclassical and mixed quantum-classical nonadiabatic methods and the reactive scattering of H2 on Ag(111) using the molecular dynamics with electronic friction method. Together, they exemplify the broad scope of the package to study effective model Hamiltonians and realistic atomistic systems.

Program Synthesis of Sparse Algorithms for Wave Function and Energy Prediction in Grid-Based Quantum Simulations.

We have recently shown how program synthesis (PS), or the concept of "self-writing code", can generate novel algorithms that solve the vibrational Schrödinger equation, providing approximations to the allowed wave functions for bound, one-dimensional (1-D) potential energy surfaces (PESs). The resulting algorithms use a grid-based representation of the underlying wave function ψ(x) and PES V(x), providing codes which represent approximations to standard discrete variable representation (DVR) methods. In this Article, we show how this inductive PS strategy can be improved and modified to enable prediction of both vibrational wave functions and energy eigenvalues of representative model PESs (both 1-D and multidimensional). We show that PS can generate algorithms that offer some improvements in energy eigenvalue accuracy over standard DVR schemes; however, we also demonstrate that PS can identify accurate numerical methods that exhibit desirable computational features, such as employing very sparse (tridiagonal) matrices. The resulting PS-generated algorithms are initially developed and tested for 1-D vibrational eigenproblems, before solution of multidimensional problems is demonstrated; we find that our new PS-generated algorithms can reduce calculation times for grid-based eigenvector computation by an order of magnitude or more. More generally, with further development and optimization, we anticipate that PS-generated algorithms based on effective Hamiltonian approximations, such as those proposed here, could be useful in direct simulations of quantum dynamics via wave function propagation and evaluation of molecular electronic structure.

Predicting Molecular Photochemistry Using Machine-Learning-Enhanced Quantum Dynamics Simulations.

The processes which occur after molecules absorb light underpin an enormous range of fundamental technologies and applications, including photocatalysis to enable new chemical transformations, sunscreens to protect against the harmful effects of UV overexposure, efficient photovoltaics for energy generation from sunlight, and fluorescent probes to image the intricate details of complex biomolecular structures. Reflecting this broad range of applications, an enormously versatile set of experiments are now regularly used to interrogate light-driven chemical dynamics, ranging from the typical ultrafast transient absorption spectroscopy used in many university laboratories to the inspiring central facilities around the world, such as the next-generation of X-ray free-electron lasers.Computer simulations of light-driven molecular and material dynamics are an essential route to analyzing the enormous amount of transient electronic and structural data produced by these experimental sources. However, to date, the direct simulation of molecular photochemistry remains a frontier challenge in computational chemical science, simultaneously demanding the accurate treatment of molecular electronic structure, nuclear dynamics, and the impact of nonadiabatic couplings.To address these important challenges and to enable new computational methods which can be integrated with state-of-the-art experimental capabilities, the past few years have seen a burst of activity in the development of "direct" quantum dynamics methods, merging the machine learning of potential energy surfaces (PESs) and nonadiabatic couplings with accurate quantum propagation schemes such as the multiconfiguration time-dependent Hartree (MCTDH) method. The result of this approach is a new generation of direct quantum dynamics tools in which PESs are generated in tandem with wave function propagation, enabling accurate "on-the-fly" simulations of molecular photochemistry. These simulations offer an alternative route toward gaining quantum dynamics insights, circumventing the challenge of generating ab initio electronic structure data for PES fitting by instead only demanding expensive energy evaluations as and when they are needed.In this Account, we describe the chronological evolution of our own contributions to this field, focusing on describing the algorithmic developments that enable direct MCTDH simulations for complex molecular systems moving on multiple coupled electronic states. Specifically, we highlight active learning strategies for generating PESs during grid-based quantum chemical dynamics simulations, and we discuss the development and impact of novel diabatization schemes to enable direct grid-based simulations of photochemical dynamics; these developments are highlighted in a series of benchmark molecular simulations of systems containing multiple nuclear degrees of freedom moving on multiple coupled electronic states. We hope that the ongoing developments reported here represent a major step forward in tools for modeling excited-state chemistry such as photodissociation, proton and electron transfer, and ultrafast energy dissipation in complex molecular systems.

Analyzing Grid-Based Direct Quantum Molecular Dynamics Using Non-Linear Dimensionality Reduction.

Grid-based schemes for simulating quantum dynamics, such as the multi-configuration time-dependent Hartree (MCTDH) method, provide highly accurate predictions of the coupled nuclear and electronic dynamics in molecular systems. Such approaches provide a multi-dimensional, time-dependent view of the system wavefunction represented on a coordinate grid; in the case of non-adiabatic simulations, additional information about the state populations adds a further layer of complexity. As such, wavepacket motion on potential energy surfaces which couple many nuclear and electronic degrees-of-freedom can be extremely challenging to analyse in order to extract physical insight beyond the usual expectation-value picture. Here, we show that non-linear dimensionality reduction (NLDR) methods, notably diffusion maps, can be adapted to extract information from grid-based wavefunction dynamics simulations, providing insight into key nuclear motions which explain the observed dynamics. This approach is demonstrated for 2-D and 9-D models of proton transfer in salicylaldimine, as well as 8-D and full 12-D simulations of cis-trans isomerization in ethene; these simulations demonstrate how NLDR can provide alternative views of wavefunction dynamics, and also highlight future developments.

Experimental and Computational Analysis of Para-Hydroxy Methylcinnamate following Photoexcitation.

Para-hydroxy methylcinnamate is part of the cinnamate family of molecules. Experimental and computational studies have suggested conflicting non-radiative decay routes after photoexcitation to its S1(ππ*) state. One non-radiative decay route involves intersystem crossing mediated by an optically dark singlet state, whilst the other involves direct intersystem crossing to a triplet state. Furthermore, irrespective of the decay mechanism, the lifetime of the initially populated S1(ππ*) state is yet to be accurately measured. In this study, we use time-resolved ion-yield and photoelectron spectroscopies to precisely determine the S1(ππ*) lifetime for the s-cis conformer of para-hydroxy methylcinnamate, combined with time-dependent density functional theory to determine the major non-radiative decay route. We find the S1(ππ*) state lifetime of s-cis para-hydroxy methylcinnamate to be ∼2.5 picoseconds, and the major non-radiative decay route to follow the [1ππ*→1nπ*→3ππ*→S0] pathway. These results also concur with previous photodynamical studies on structurally similar molecules, such as para-coumaric acid and methylcinnamate.

Solving the Schrödinger equation using program synthesis.

We demonstrate that a program synthesis approach based on a linear code representation can be used to generate algorithms that approximate the ground-state solutions of one-dimensional time-independent Schrödinger equations constructed with bound polynomial potential energy surfaces (PESs). Here, an algorithm is constructed as a linear series of instructions operating on a set of input vectors, matrices, and constants that define the problem characteristics, such as the PES. Discrete optimization is performed using simulated annealing in order to identify sequences of code-lines, operating on the program inputs that can reproduce the expected ground-state wavefunctions ψ(x) for a set of target PESs. The outcome of this optimization is not simply a mathematical function approximating ψ(x) but is, instead, a complete algorithm that converts the input vectors describing the system into a ground-state solution of the Schrödinger equation. These initial results point the way toward an alternative route for developing novel algorithms for quantum chemistry applications.

Identifying Barrierless Mechanisms for Benzene Formation in the Interstellar Medium Using Permutationally Invariant Reaction Discovery.

Complex chemical reaction environments, such as those found in combustion engines, the upper atmosphere, or the interstellar medium, can contain large numbers of different reactive species participating in similarly large numbers of different chemical reactions. In such settings, identifying the most-likely multistep reaction mechanisms which lead to the production of a particular defined product species is an extremely challenging problem, requiring search and evaluation over a large number of different possible candidate mechanisms while also addressing the permutational challenges posed when considering a large number of reaction routes available to sets of identical molecular species. In this article, the problem of generating candidate reaction mechanisms which form a defined product from a diverse set of reactive molecules is cast as a discrete optimization of a permutationally invariant cost function describing similarity between the target product and the product generated by a trial reaction mechanism. This approach is demonstrated by generating 2230 candidate reaction mechanisms which form benzene from diverse sets of reactive molecules which have been experimentally identified in the interstellar medium. By screening this set of autogenerated mechanisms, using dispersion-corrected DFT to evaluate reaction energies and activation barriers, we identify several candidate barrierless reaction mechanisms (both previously proposed and new) for benzene formation which may operate in the low temperatures found in the interstellar medium and could be investigated further to supplement existing microkinetic models.

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.

Direct Grid-Based Nonadiabatic Dynamics on Machine-Learned Potential Energy Surfaces: Application to Spin-Forbidden Processes.

We have recently shown how high-accuracy wave function grid-based propagation schemes, such as the multiconfiguration time-dependent Hartree (MCTDH) method, can be combined with machine-learning (ML) descriptions of PESs to yield an "on-the-fly" direct dynamics scheme which circumvents potential energy surface (PES) prefitting. To date, our approach has been demonstrated in the ground-state dynamics and nonadiabatic spin-allowed dynamics of several molecular systems. Expanding on this successful previous work, this Article demonstrates how our ML-based quantum dynamics scheme can be adapted to model nonadiabatic dynamics for spin-forbidden processes such as intersystem crossing (ISC), opening up new possibilities for modeling chemical dynamic phenomena driven by spin-orbit coupling. After describing modifications to diabatization schemes to enable accurate and robust treatment or electronic states of different spin-multiplicity, we demonstrate our methodology in applications to modeling ISC in SO2 and thioformaldehyde, benchmarking our results against previous trajectory- and grid-based calculations. As a relatively efficient tool for modeling spin-forbidden nonadiabatic dynamics without demanding any prefitting of PESs, our overall strategy is a potentially powerful tool for modeling important photochemical systems, such as photoactivated pro-drugs and organometallic catalysts.

A new diabatization scheme for direct quantum dynamics: Procrustes diabatization.

We present a new scheme for diabatizing electronic potential energy surfaces for use within the recently implemented direct-dynamics grid-based class of computational nuclear quantum dynamics methods, called Procrustes diabatization. Calculations on the well-studied molecular systems LiF and the butatriene cation, using both Procrustes diabatization and the previously implemented propagation and projection diabatization schemes, have allowed detailed comparisons to be made, which indicate that the new method combines the best features of the older approaches; it generates smooth surfaces, which cross at the correct molecular geometries, reproduces interstate couplings accurately, and hence allows the correct modeling of non-adiabatic dynamics.