A size-consistent multi-state mapping approach to surface hopping.

We develop a multi-state generalization of the recently proposed mapping approach to surface hopping (MASH) for the simulation of electronically nonadiabatic dynamics. This new approach extends the original MASH method to be able to treat systems with more than two electronic states. It differs from previous approaches in that it is size consistent and rigorously recovers the original two-state MASH in the appropriate limits. We demonstrate the accuracy of the method by applying it to a series of model systems for which exact benchmark results are available, and we find that the method is well suited to the simulation of photochemical relaxation processes.

A MASH simulation of the photoexcited dynamics of cyclobutanone.

In response to a community prediction challenge, we simulate the nonadiabatic dynamics of cyclobutanone using the mapping approach to surface hopping (MASH). We consider the first 500 fs of relaxation following photoexcitation to the S2 state and predict the corresponding time-resolved electron-diffraction signal that will be measured by the planned experiment. 397 ab initio trajectories were obtained on the fly with state-averaged complete active space self-consistent field using a (12,11) active space. To obtain an estimate of the potential systematic error, 198 of the trajectories were calculated using an aug-cc-pVDZ basis set and 199 with a 6-31+G* basis set. MASH is a recently proposed independent trajectory method for simulating nonadiabatic dynamics, originally derived for two-state problems. As there are three relevant electronic states in this system, we used a newly developed multi-state generalization of MASH for the simulation: the uncoupled spheres multi-state MASH method (unSMASH). This study, therefore, serves both as an investigation of the photodissociation dynamics of cyclobutanone, and also as a demonstration of the applicability of unSMASH to ab initio simulations. In line with previous experimental studies, we observe that the simulated dynamics is dominated by three sets of dissociation products, C3H6 + CO, C2H4 + C2H2O, and C2H4 + CH2 + CO, and we interpret our predicted electron-diffraction signal in terms of the key features of the associated dissociation pathways.

Perturbatively corrected ring-polymer instanton theory for accurate tunneling splittings.

We introduce an approach for calculating perturbative corrections to the ring-polymer instanton approximation to tunneling splittings (RPI+PC) by computing higher-order terms in the asymptotic expansion in ℏ. The resulting method goes beyond standard instanton theory by using information on the third and fourth derivatives of the potential along the tunneling path to include additional anharmonic effects. This leads to significant improvements both in systems with low barriers and in systems with anharmonic modes. We demonstrate the applicability of RPI+PC to molecular systems by computing the tunneling splitting in full-dimensional malonaldehyde and a deuterated derivative. Comparing to both experiment and recent quantum mechanical benchmark results, we find that our perturbative correction reduces the error from -11% to 2% for hydrogen transfer and performs even better for the deuterated case. This makes our approach more accurate than previous calculations using diffusion Monte Carlo and path-integral molecular dynamics while being more computationally efficient.

Vibrational strong coupling in liquid water from cavity molecular dynamics.

We assess the cavity molecular dynamics method for the calculation of vibrational polariton spectra using liquid water as a specific example. We begin by disputing a recent suggestion that nuclear quantum effects may lead to a broadening of polariton bands, finding instead that they merely result in anharmonic red shifts in the polariton frequencies. We go on to show that our simulated cavity spectra can be reproduced to graphical accuracy with a harmonic model that uses just the cavity-free spectrum and the geometry of the cavity as input. We end by showing that this harmonic model can be combined with the experimental cavity-free spectrum to give results in good agreement with optical cavity measurements. Since the input to our harmonic model is equivalent to the input to the transfer matrix method of applied optics, we conclude that cavity molecular dynamics cannot provide any more insight into the effect of vibrational strong coupling on the absorption spectrum than this transfer matrix method, which is already widely used by experimentalists to corroborate their cavity results.

Improved microcanonical instanton theory.

Canonical (thermal) instanton theory is now routinely applicable to complex gas-phase reactions and allows for the accurate description of tunnelling in highly non-separable systems. Microcanonical instanton theory is by contrast far less well established. Here, we demonstrate that the best established microcanonical theory [S. Chapman, B. C. Garrett and W. H. Miller, J. Chem. Phys., 1975, 63, 2710-2716], fails to accurately describe the deep-tunnelling regime for systems where the frequencies of the orthogonal modes change rapidly along the instanton path. By taking a first principles approach to the derivation of microcanonical instanton theory, we obtain an improved method, which accurately recovers the thermal instanton rate when integrated over energy. The resulting theory also correctly recovers the separable limit and can be thought of as an instanton generalisation of Rice-Ramsperger-Kassel-Marcus (RRKM) theory. When combined with the density-of-states approach [W. Fang, P. Winter and J. O. Richardson, J. Chem. Theory Comput., 2021, 17, 40-55], this new method can be straightforwardly applied to real molecular systems.

Fast quasi-centroid molecular dynamics.

We describe a fast implementation of the quasi-centroid molecular dynamics (QCMD) method in which the quasi-centroid potential of mean force is approximated as a separable correction to the classical interaction potential. This correction is obtained by first calculating quasi-centroid radial and angular distribution functions in a short path integral molecular dynamics simulation and then using iterative Boltzmann inversion to obtain an effective classical potential that reproduces these distribution functions in a classical NVT simulation. We illustrate this approach with example applications to the vibrational spectra of gas phase molecules, obtaining excellent agreement with QCMD reference calculations for water and ammonia and good agreement with the quantum mechanical vibrational spectrum of methane.

Confirming the role of nuclear tunneling in aqueous ferrous-ferric electron transfer.

We revisit the well-known aqueous ferrous-ferric electron transfer reaction in order to address recent suggestions that nuclear tunneling can lead to significant deviation from the linear response assumption inherent in the Marcus picture of electron transfer. A recent study of this reaction by Richardson and co-workers [Phys. Chem. Chem. Phys. 22, 10687 (2020)] has found a large difference between their new path-integral method, golden-rule quantum transition state theory (GR-QTST), and the saddle point approximation of Wolynes (Wolynes theory). They suggested that this difference could be attributed to the existence of multiple tunneling pathways, leading Wolynes theory to significantly overestimate the rate. This was used to argue that the linear response assumptions of Marcus theory may break down for liquid systems when tunneling is important. If true, this would imply that the commonly used method for studying such systems, where the problem is mapped onto a spin-boson model, is invalid. However, we have recently shown that size inconsistency in GR-QTST can lead to poor predictions of the rate in systems with many degrees of freedom. We have also suggested an improved method, the path-integral linear golden-rule (LGR) approximation, which fixes this problem. Here, we demonstrate that the GR-QTST results for ferrous-ferric electron transfer are indeed dominated by its size consistency error. Furthermore, by comparing the LGR and Wolynes theory results, we confirm the established picture of nuclear tunneling in this system. Finally, by comparing our path-integral results to those obtained by mapping onto the spin-boson model, we reassess the importance of anharmonic effects and the accuracy of this commonly used mapping approach.

An improved path-integral method for golden-rule rates.

We present a simple method for the calculation of reaction rates in the Fermi golden-rule limit, which accurately captures the effects of tunneling and zero-point energy. The method is based on a modification of the recently proposed golden-rule quantum transition state theory (GR-QTST) of Thapa, Fang, and Richardson [J. Chem. Phys. 150, 104107 (2019)]. While GR-QTST is not size consistent, leading to the possibility of unbounded errors in the rate, our modified method has no such issue and so can be reliably applied to condensed phase systems. Both methods involve path-integral sampling in a constrained ensemble; the two methods differ, however, in the choice of constraint functional. We demonstrate numerically that our modified method is as accurate as GR-QTST for the one-dimensional model considered by Thapa and co-workers. We then study a multidimensional spin-boson model, for which our method accurately predicts the true quantum rate, while GR-QTST breaks down with an increasing number of boson modes in the discretization of the spectral density. Our method is able to accurately predict reaction rates in the Marcus inverted regime without the need for the analytic continuation required by Wolynes theory.

A general non-adiabatic quantum instanton approximation.

We present a general quantum instanton approach to calculating reaction rates for systems with two electronic states and arbitrary values of the electronic coupling. This new approach, which we call the non-adiabatic quantum instanton (NAQI) approximation, reduces to Wolynes theory in the golden rule limit and to a recently proposed projected quantum instanton method in the adiabatic limit. As in both of these earlier theories, the NAQI approach is based on making a saddle point approximation to the time integral of a reactive flux autocorrelation function, although with a generalized definition of the projection operator onto the product states. We illustrate the accuracy of the approach by comparison with exact rates for one dimensional scattering problems and discuss its applicability to more complex reactions.

Path integral methods for reaction rates in complex systems.

We shall use this introduction to the Faraday Discussion on quantum effects in complex systems to review the recent progress that has been made in using imaginary time path integral methods to calculate chemical reaction rates. As a result of this progress, it is now routinely possible to calculate accurate rate constants including quantum mechanical zero point energy and tunnelling effects for arbitrarily complex (anharmonic and multi-dimensional) systems. This is true in the adiabatic (Born-Oppenheimer) limit, in the non-adiabatic (Fermi Golden Rule) limit, and everywhere between these two limits in the normal Marcus regime. Quantum mechanical effects on reaction rates can be enormous, even at room temperature, and the problem of including these effects in simulations of a wide variety of chemical reactions in complex systems has now effectively been solved.

An analysis of isomorphic RPMD in the golden rule limit.

We analyze the golden rule limit of the recently proposed isomorphic ring polymer (iso-RP) method. This method aims to combine an exact expression for the quantum mechanical partition function of a system with multiple electronic states with a pre-existing mixed quantum-classical (MQC) dynamics approximation, such as fewest switches surface hopping. Since the choice of the MQC method adds a degree of flexibility, we simplify the analysis by assuming that the dynamics used correctly reproduces the exact golden rule rate for a nonadiabatic (e.g., electron transfer) reaction in the high temperature limit. Having made this assumption, we obtain an expression for the iso-RP rate in the golden rule limit that is valid at any temperature. We then compare this rate with the exact rate for a series of simple spin-boson models. We find that the iso-RP method does not correctly predict how nuclear quantum effects affect the reaction rate in the golden rule limit. Most notably, it does not capture the quantum asymmetry in a conventional (Marcus) plot of the logarithm of the reaction rate against the thermodynamic driving force, and it also significantly overestimates the correct quantum mechanical golden rule rate for activationless electron transfer reactions. These results are analyzed and their implications discussed for the applicability of the iso-RP method to more general nonadiabatic reactions.

On the calculation of quantum mechanical electron transfer rates.

We present a simple interpolation formula for the rate of an electron transfer reaction as a function of the electronic coupling strength. The formula only requires the calculation of Fermi golden rule and Born-Oppenheimer rates and so can be combined with any methods that are able to calculate these rates. We first demonstrate the accuracy of the formula by applying it to a one dimensional scattering problem for which the exact quantum mechanical, Fermi golden rule, and Born-Oppenheimer rates are readily calculated. We then describe how the formula can be combined with the Wolynes theory approximation to the golden rule rate, and the ring polymer molecular dynamics (RPMD) approximation to the Born-Oppenheimer rate, and used to capture the effects of nuclear tunneling, zero point energy, and solvent friction on condensed phase electron transfer reactions. Comparison with exact hierarchical equations of motion results for a demanding set of spin-boson models shows that the interpolation formula has an error comparable to that of RPMD rate theory in the adiabatic limit, and that of Wolynes theory in the nonadiabatic limit, and is therefore as accurate as any method could possibly be that attempts to generalize these methods to arbitrary electronic coupling strengths.

Analytic continuation of Wolynes theory into the Marcus inverted regime.

The Wolynes theory of electronically nonadiabatic reaction rates [P. G. Wolynes, J. Chem. Phys. 87, 6559 (1987)] is based on a saddle point approximation to the time integral of a reactive flux autocorrelation function in the nonadiabatic (golden rule) limit. The dominant saddle point is on the imaginary time axis at tsp=iλspℏ, and provided λsp lies in the range -ß/2≤λsp≤ß/2, it is straightforward to evaluate the rate constant using information obtained from an imaginary time path integral calculation. However, if λsp lies outside this range, as it does in the Marcus inverted regime, the path integral diverges. This has led to claims in the literature that Wolynes theory cannot describe the correct behaviour in the inverted regime. Here we show how the imaginary time correlation function obtained from a path integral calculation can be analytically continued to λsp<-ß/2, and the continuation used to evaluate the rate in the inverted regime. Comparison with exact golden rule results for a spin-boson model and a more demanding (asymmetric and anharmonic) model of electronic predissociation shows that the theory is just as accurate in the inverted regime as it is in the normal regime.

Magnetoelectroluminescence in organic light-emitting diodes.

The magnetoelectroluminescence of conjugated organic polymer films is widely accepted to arise from a polaron pair mechanism, but their magnetoconductance is less well understood. Here we derive a new relationship between the experimentally measurable magnetoelectroluminescence and magnetoconductance and the theoretically calculable singlet yield of the polaron pair recombination reaction. This relationship is expected to be valid regardless of the mechanism of the magnetoconductance, provided the mobilities of the free polarons are independent of the applied magnetic field (i.e., provided one discounts the possibility of spin-dependent transport). We also discuss the semiclassical calculation of the singlet yield of the polaron pair recombination reaction for materials such as poly(2,5-dioctyloxy-paraphenylene vinylene) (DOO-PPV), the hyperfine fields in the polarons of which can be extracted from light-induced electron spin resonance measurements. The resulting theory is shown to give good agreement with experimental data for both normal (H-) and deuterated (D-) DOO-PPV over a wide range of magnetic field strengths once singlet-triplet dephasing is taken into account. Without this effect, which has not been included in any previous simulation of magnetoelectroluminescence, it is not possible to reproduce the experimental data for both isotopologues in a consistent fashion. Our results also indicate that the magnetoconductance of DOO-PPV cannot be solely due to the effect of the magnetic field on the dissociation of polaron pairs.

Recovering Marcus Theory Rates and Beyond without the Need for Decoherence Corrections: The Mapping Approach to Surface Hopping.

It is well-known that fewest-switches surface hopping (FSSH) fails to correctly capture the quadratic scaling of rate constants with diabatic coupling in the weak-coupling limit, as expected from Fermi's golden rule and Marcus theory. To address this deficiency, the most widely used approach is to introduce a "decoherence correction", which removes the inconsistency between the wave function coefficients and the active state. Here we investigate the behavior of a new nonadiabatic trajectory method, called the mapping approach to surface hopping (MASH), on systems that exhibit an incoherent rate behavior. Unlike FSSH, MASH hops between active surfaces deterministically and can never have an inconsistency between the wave function coefficients and the active state. We show that MASH not only can describe rates for intermediate and strong diabatic coupling but also can accurately reproduce the results of Marcus theory in the golden-rule limit, without the need for a decoherence correction. MASH is therefore a significant improvement over FSSH in the simulation of nonadiabatic reactions.

How Quantum is the Resonance Behavior in Vibrational Polariton Chemistry?

Recent experiments in polariton chemistry have demonstrated that reaction rates can be modified by vibrational strong coupling to an optical cavity mode. Importantly, this modification occurs only when the frequency of the cavity mode is tuned to closely match a molecular vibrational frequency. This sharp resonance behavior has proved to be difficult to capture theoretically. Only recently did Lindoy et al. [ Nat. Commun. 2023, 14, 2733] report the first instance of a sharp resonant effect in the cavity-modified rate simulated in a model system using exact quantum dynamics. We investigate the same model system with a different method, ring-polymer molecular dynamics (RPMD), which captures quantum statistics but treats dynamics classically. We find that RPMD does not reproduce this sharp resonant feature at the well frequency, and we discuss the implications of this finding for future studies of vibrational polariton chemistry.

Fast Quasi-Centroid Molecular Dynamics for Water and Ice.

We describe how the fast quasi-centroid molecular dynamics (f-QCMD) method can be applied to condensed-phase systems by approximating the quasi-centroid potential of mean force as a sum of inter- and intramolecular corrections to the classical interaction potential. The corrections are found by using a regularized iterative Boltzmann inversion procedure to recover the inter- and intramolecular quasi-centroid distribution functions obtained from a path integral molecular dynamics simulation. The resulting methodology is found to give good agreement with a previously published QCMD dipole absorption spectrum for liquid water and satisfactory agreement for ice. It also gives good agreement with spectra from a recent implementation of CMD that uses a precomputed elevated temperature potential of mean force. Modern centroid molecular dynamics methods, therefore, appear to be reaching a consensus regarding the impact of nuclear quantum effects on the vibrational spectra of water and ice.