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.

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.

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.

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.

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.

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.

Can we use on-the-fly quantum simulations to connect molecular structure and sunscreen action?

The ground-up design of new molecular sunscreens, with improved photostability, absorbance and spectral coverage, stands as a challenge to fundamental chemical science. Correlating sunscreen molecular structure and function requires detailed insight into the relaxation pathways available following photoexcitation; however, the complex coupled electron/nuclear dynamics in these systems stands as a tough challenge to computational chemistry. To address this challenge, we have recently developed efficient and accurate simulation methods to model non-adiabatic dynamics of general molecular systems as a route to correlating photoinduced dynamics and potential sunscreen activity. Our approach, combining the multi-configuration time-dependent Hartree (MCTDH) method with PESs generated using machine-learning, represents a new "on-the-fly" strategy for accurate wavefunction propagation, with the potential to provide a new "black box" strategy for interrogating ultrafast dynamics in general photoinduced energy transport processes. Here, we illustrate our attempts to apply this methodology to study the ultrafast photochemistry of four compounds derived from mycosporine-like amino acids (MAAs), compounds which are believed to act as microbial sunscreens in micro-organisms such as algae. Specifically, we investigate how the choice of active vibrational space and diabatic electronic states in MCTDH strongly influences the predictions of ultrafast dynamic relaxation in representative MAAs. Our results serve to demonstrate that "on-the-fly" quantum dynamics using MCTDH is increasingly viable, but important barriers relating to coordinate choices and diabatisation remain.

Automatic Proposal of Multistep Reaction Mechanisms using a Graph-Driven Search.

Proposing and testing mechanistic hypotheses stands as one of the key applications of contemporary computational chemistry. In the majority of computational mechanistic analyses, the individual elementary steps leading from reactants to products are proposed by the user, based on learned chemical knowledge, intuition, or comparison to an existing well-characterized mechanism for a closely related chemical reaction. However, the prerequisite of prior chemical knowledge is a barrier to automated (or "black box") mechanistic generation and assessment, and it may simultaneously preclude mechanistic proposals that lie outside the "standard" chemical reaction set. In this Article, we propose a simple random-walk algorithm that searches for the set of elementary chemical reactions that transform defined reactant structures into target products. Our approach operates exclusively in the space of molecular connectivity matrices, seeking the set of chemically sensible bonding changes that link connectivity matrices for input reactant and product structures. We subsequently illustrate how atomic coordinates for each elementary reaction can be generated under the action of a graph-restraining potential, prior to further analysis by quantum chemical calculations. Our approach is successfully demonstrated for carbon monoxide oxidation, the water-gas shift reaction, and n-hexane aromatization, all catalyzed by Pt nanoparticles.

The role of nuclear quantum effects in the relative stability of hexagonal and cubic ice.

At atmospheric pressure, hexagonal ice (Ih) is thermodynamically stable relative to cubic ice (Ic), although the magnitude and underlying physical origin of this stability difference are not well defined. Pure Ic crystals are not accessible experimentally, and hence computer simulations have often been used to interrogate the relative stabilities of Ih and Ic; however, these simulations are dominated by molecular interaction models that ignore the intramolecular flexibility of individual water molecules, do not describe intermolecular hydrogen-bonding with sufficient accuracy, or ignore the role of nuclear quantum effects (NQEs) such as zero-point energy. Here, we show that when comparing the relative stability of Ih and Ic using a flexible, anharmonic molecular interaction model, while also accurately accounting for NQEs, a new picture emerges: Ih is stabilized relative to Ic as a result of subtle differences in the intramolecular geometries and intermolecular interactions of water molecules which are modulated by NQEs. Our simulations hence suggest that NQEs are a major contributor to the stabilization of Ih under terrestrial conditions and thus contribute to the well-known hexagonal (sixfold) symmetry of ice crystals.

Direct quantum dynamics using variational Gaussian wavepackets and Gaussian process regression.

The method of direct variational quantum nuclear dynamics in a basis of Gaussian wavepackets, combined with the potential energy surfaces fitted on-the-fly using Gaussian process regression, is described together with its implementation. Enabling exact and efficient analytic evaluation of Hamiltonian matrix elements, this approach allows for black-box quantum dynamics of multidimensional anharmonic molecular systems. Example calculations of intra-molecular proton transfer on the electronic ground state of salicylaldimine are provided, and future algorithmic improvements as well as the potential for multiple-state non-adiabatic dynamics are discussed.

MCTDH on-the-fly: Efficient grid-based quantum dynamics without pre-computed potential energy surfaces.

We present significant algorithmic improvements to a recently proposed direct quantum dynamics method, based upon combining well established grid-based quantum dynamics approaches and expansions of the potential energy operator in terms of a weighted sum of Gaussian functions. Specifically, using a sum of low-dimensional Gaussian functions to represent the potential energy surface (PES), combined with a secondary fitting of the PES using singular value decomposition, we show how standard grid-based quantum dynamics methods can be dramatically accelerated without loss of accuracy. This is demonstrated by on-the-fly simulations (using both standard grid-based methods and multi-configuration time-dependent Hartree) of both proton transfer on the electronic ground state of salicylaldimine and the non-adiabatic dynamics of pyrazine.

Harmonic-phase path-integral approximation of thermal quantum correlation functions.

We present an approximation to the thermal symmetric form of the quantum time-correlation function in the standard position path-integral representation. By transforming to a sum-and-difference position representation and then Taylor-expanding the potential energy surface of the system to second order, the resulting expression provides a harmonic weighting function that approximately recovers the contribution of the phase to the time-correlation function. This method is readily implemented in a Monte Carlo sampling scheme and provides exact results for harmonic potentials (for both linear and non-linear operators) and near-quantitative results for anharmonic systems for low temperatures and times that are likely to be relevant to condensed phase experiments. This article focuses on one-dimensional examples to provide insights into convergence and sampling properties, and we also discuss how this approximation method may be extended to many-dimensional systems.

Photosynthesis, pigment-protein complexes and electronic energy transport: simple models for complicated processes.

In this review, we discuss our recent work on modelling biological pigment-protein complexes, such as the Fenna-Matthews-Olson complex and light-harvesting complex-II, to explain their electronic energy transport properties. In particular, we highlight how a network-based analysis approach, where the light-absorbing pigments are treated as a network of interconnected nodes, can provide a qualitative picture of quantum dynamic energy transport. With this in mind, we demonstrate how other properties such as robustness to environmental changes can be assessed in a simple and computationally tractable manner. Such analyses could prove useful for the design of artificial energy transport networks such as those which might find application in solar cells.

Accelerated path-integral simulations using ring-polymer interpolation.

Imaginary-time path-integral (PI) molecular simulations can be used to calculate exact quantum statistical mechanical properties for complex systems containing many interacting atoms and molecules. The limiting computational factor in a PI simulation is typically the evaluation of the potential energy surface (PES) and forces at each ring-polymer "bead"; for an n-bead ring-polymer, a PI simulation is typically n times greater than the corresponding classical simulation. To address the increased computational effort of PI simulations, several approaches have been developed recently, most notably based on the idea of ring-polymer contraction which exploits either the separation of the PES into short-range and long-range contributions or the availability of a computationally inexpensive PES which can be incorporated to effectively smooth the ring-polymer PES; neither approach is satisfactory in applications to systems modeled by PESs given by on-the-fly ab initio calculations. In this article, we describe a new method, ring-polymer interpolation (RPI), which can be used to accelerate PI simulations without any prior assumptions about the PES. In simulations of liquid water modeled by an empirical PES (or force field) under ambient conditions, where quantum effects are known to play a subtle role in influencing experimental observables such as radial distribution functions, we find that RPI can accurately reproduce the results of fully-converged PI simulations, albeit with far fewer PES evaluations. This approach therefore opens the possibility of large-scale PI simulations using ab initio PESs evaluated on-the-fly without the drawbacks of current methods.

Efficient and accurate evaluation of potential energy matrix elements for quantum dynamics using Gaussian process regression.

Solution of the time-dependent Schrödinger equation using a linear combination of basis functions, such as Gaussian wavepackets (GWPs), requires costly evaluation of integrals over the entire potential energy surface (PES) of the system. The standard approach, motivated by computational tractability for direct dynamics, is to approximate the PES with a second order Taylor expansion, for example centred at each GWP. In this article, we propose an alternative method for approximating PES matrix elements based on PES interpolation using Gaussian process regression (GPR). Our GPR scheme requires only single-point evaluations of the PES at a limited number of configurations in each time-step; the necessity of performing often-expensive evaluations of the Hessian matrix is completely avoided. In applications to 2-, 5-, and 10-dimensional benchmark models describing a tunnelling coordinate coupled non-linearly to a set of harmonic oscillators, we find that our GPR method results in PES matrix elements for which the average error is, in the best case, two orders-of-magnitude smaller and, in the worst case, directly comparable to that determined by any other Taylor expansion method, without requiring additional PES evaluations or Hessian matrices. Given the computational simplicity of GPR, as well as the opportunities for further refinement of the procedure highlighted herein, we argue that our GPR methodology should replace methods for evaluating PES matrix elements using Taylor expansions in quantum dynamics simulations.