Branching corrected surface hopping: Resetting wavefunction coefficients based on judgement of wave packet reflection

Detailed Complementary Consistency: Wave Function Tells Particle How to Hop, Particle Tells Wave Function How to Collapse

In mixed quantum-classical dynamics, the quantum subsystem can have both wave function and particle-like descriptions. However, they may yield inconsistent results for the expectation value of the same physical quantity. We here propose a novel detailed complementary consistency (DCC) method based on the principle of detailed internal consistency. Namely, the wave function along each trajectory tells the particle how to hop, while the particle tells the wave function how to collapse based on active states in the trajectory ensemble. As benchmarked in a diverse array of representative models with localized nonadiabatic couplings, DCC not only achieves fully consistent results (i.e., identical populations calculated based on wave functions and active states) but also closely reproduces the exact quantum results. Due to the high performance, our new DCC method has great potential to give a consistent and accurate mixed quantum-classical description of general nonadiabatic dynamics after further development.

Nonadiabatic Dynamics with Exact Factorization: Implementation and Assessment

In this work, we report our implementation of several independent-trajectory mixed-quantum-classical (ITMQC) nonadiabatic dynamics methods based on exact factorization (XF) in the Libra package for nonadiabatic and excited-state dynamics. Namely, the exact factorization surface hopping (SHXF), mixed quantum-classical dynamics (MQCXF), and mean-field (MFXF) are introduced. Performance of these methods is compared to that of several traditional surface hopping schemes, such as the fewest-switches surface hopping (FSSH), branching-corrected surface hopping (BCSH), and the simplified decay of mixing (SDM), as well as conventional Ehrenfest (mean-field, MF) method. Based on a comprehensive set of 1D model Hamiltonians, we find the ranking SHXF ≈ MQCXF > BCSH > SDM > FSSH ≫ MF, with the BCSH sometimes outperforming the XF methods in terms of describing coherences. Although the MFXF method can yield reasonable populations and coherences for some cases, it does not conserve the total energy and is therefore not recommended. We also find that the branching correction for auxiliary trajectories is important for the XF methods to yield accurate populations and coherences. However, the branching correction can worsen the quality of the energy conservation in the MQCXF. Finally, we find that using the time-dependent Gaussian width approximation used in the XF methods for computing decoherence correction can improve the quality of energy conservation in the MQCXF dynamics. The parameter-free scheme of Subotnik for computing the Gaussian widths is found to deliver the best performance in situations where such widths are not known a priori.

Consistent Construction of the Density Matrix from Surface Hopping Trajectories

Proper construction of the density matrix based on surface hopping trajectories remains a difficult problem. Due to the well-known overcoherence in traditional surface hopping simulations, the electronic wave function cannot be used directly. In this work, we propose a consistent density matrix construction method, which takes the advantage of occupation of active states to rescale the coherence calculated by wave functions and ensures the intrinsic consistency of the density matrix. This new trajectory analysis method can be used for both Tully's fewest switches surface hopping (FSSH) and our recently proposed branching corrected surface hopping (BCSH). As benchmarked in both one- and two-dimensional standard scattering models, the new approach combined with BCSH trajectories achieves highly accurate time-dependent spatial distributions of adiabatic populations and coherence compared to exact quantum results.

Surface Hopping with Reliable Wave Function by Introducing Auxiliary Wave Packets to Trajectory Branching

It is well-known that the widely utilized fewest switches surface hopping method suffers from the severe overcoherence problem, and thus adiabatic populations calculated by wave functions are generally inferior to those based on active states. More importantly, to achieve a complete description of nonadiabatic dynamics, the density matrix is essential. In this paper, we present an auxiliary branching corrected surface hopping (A-BCSH) method that introduces auxiliary wave packets (WPs) on the adiabatic potential energy surfaces for trajectory branching. Both rapid and gradual separation of WP components on different surfaces are characterized, and thus the correct decoherence time along each trajectory is captured. As demonstrated in the three standard Tully models, A-BCSH exhibits excellent internal consistency. Namely, close adiabatic populations are obtained based on both wave functions and active states. In particular, A-BCSH successfully obtains a reliable time-dependent spatial distribution of the density matrix, which relies only on electronic wave functions. Due to its high performance, our A-BCSH method provides a new and highly promising perspective on further development of more consistent surface hopping with reliable wave function.

Investigating the Photodynamics of trans-Azobenzene with Coupled Trajectories

In this work, we present the first implementation of coupled-trajectory Tully surface hopping (CT-TSH) suitable for applications to molecular systems. We combine CT-TSH with the semiempirical floating occupation molecular orbital-configuration interaction electronic structure method to investigate the photoisomerization dynamics of trans-azobenzene. Our study shows that CT-TSH can capture correctly decoherence effects in this system, yielding consistent electronic and nuclear dynamics in agreement with (standard) decoherence-corrected TSH. Specifically, CT-TSH is derived from the exact factorization and the electronic coefficients' evolution is directly influenced by the coupling of trajectories, resulting in the improvement of internal consistency if compared to standard TSH.

Energy-Conserving and Thermally Corrected Neglect of Back-Reaction Approximation Method for Nonadiabatic Molecular Dynamics

In this work, the energy-conserving and thermally corrected neglect of the back-reaction approximation approach for nonadiabatic molecular dynamics in extended atomistic systems is developed. The new approach introduces three key corrections to the original method: (1) it enforces the total energy conservation, (2) it introduces an explicit coupling of the system to its environment, and (3) it introduces a renormalization of nonadiabatic couplings to account for a difference between the instantaneous nuclear kinetic energy and the kinetic energy of guiding trajectories. In the new approach, an auxiliary kinetic energy variable is introduced as an independent dynamical variable. The new approach produces nonzero equilibrium populations, whereas the original neglect of the back-reaction approximation method does not. It yields population relaxation time scales that are favorably comparable to the reference values, and it introduces an explicit and controllable way of dissipating energy into a bath without an assumption of the bath being at equilibrium.

ChemPotPy: A Python Library for Analytic Representations of Potential Energy Surfaces and Diabatic Potential Energy Matrices

Constructing analytic representations of global and semiglobal potential energy surfaces is difficult and can be laborious, and it is even harder when one needs coupled potential energy surfaces and their electronically nonadiabatic couplings. When accomplished, however, the resulting potential functions are a valuable resource. To facilitate the convenient use of potentials that have been developed, we provide a collection of existing surfaces in a library with consistent units and formats. A potential energy surface library of this type, namely PotLib, was built more than 20 years ago. However, that library only provided pristine Fortran subroutines for each potential energy surface, and therefore, it is not as user-friendly as would be desirable. Here, we report the creation of ChemPotPy, a CHEMical library of POTential energy surfaces in PYthon. ChemPotPy is a user-friendly library for analytic representation of single-state and multistate potential energy surfaces and couplings. A given entry in the library contains an analytic potential energy function or analytic functions for a set of coupled potential energy surfaces, and depending on the case, it may also include analytic or numerical gradients, nonadiabatic coupling vectors, and/or diabatic potential energy matrices and their gradients. Only three inputs, namely, the chemical formula of the system, the name of the potential energy surface or surface set, and the Cartesian geometry, are required. ChemPotPy uses the same units for input and output quantities of all surfaces and surface sets to facilitate general interfaces with the dynamics programs. The initial version of the library contains 338 entries, and we anticipate that more will be added in the future.

What is Missing in the Mean Field Description of Spatial Distribution of Population? Important Role of Auxiliary Wave Packets in Trajectory Branching

When the traditional Ehrenfest mean field approach is employed to simulate nonadiabatic dynamics, an effective wave packet (WP) on the average potential energy surface (PES) is utilized to describe the nuclear motion. In the fully quantum picture, however, the WP components on different adiabatic PESs gradually separate in space because they evolve under different velocities and forces. Due to trajectory branching of the WP components, proper decoherence needs to be taken into account, and the spatial distribution of population cannot be described by a single effective WP. Here, we propose an auxiliary branching corrected mean field (A-BCMF) method, where trajectories of auxiliary WPs on adiabatic PESs are introduced. As benchmarked in the three standard Tully models, A-BCMF not only gives correct channel populations but also captures an accurate time-dependent spatial distribution of population. Thereby, we reveal the important role of auxiliary WPs in solving intrinsic problems of the widely used mean field description of nonadiabatic dynamics.

Different flavors of exact-factorization-based mixed quantum-classical methods for multistate dynamics

The exact factorization approach has led to the development of new mixed quantum-classical methods for simulating coupled electron-ion dynamics. We compare their performance for dynamics when more than two electronic states are occupied at a given time, and analyze: (1) the use of coupled versus auxiliary trajectories in evaluating the electron-nuclear correlation terms, (2) the approximation of using these terms within surface-hopping and Ehrenfest frameworks, and (3) the relevance of the exact conditions of zero population transfer away from nonadiabatic coupling regions and total energy conservation. Dynamics through the three-state conical intersection in the uracil radical cation as well as polaritonic models in one dimension are studied.

Learning Decoherence Time Formulas for Surface Hopping from Quantum Dynamics

Surface hopping simulations have achieved great success in many different fields, but their reliability has long been limited by the overcoherence problem. We here present a general machine learning assisted approach to identify optimal decoherence time formulas for surface hopping using exact quantum dynamics as references. In order to avoid computationally expensive force calculations, we use the nuclear kinetic energy and the adiabatic energy difference to iteratively generate the descriptor space. Through multilayer screening of the candidate descriptors and discrete optimization of the relevant parameters, we obtain new energy-based decoherence time formulas. As benchmarked in thousands of diverse multilevel systems and six standard scattering models, surface hopping with our new decoherence time formulas nearly reproduces the exact quantum dynamics while maintaining high efficiency. Thereby, our approach provides a promising avenue for systematically improving the accuracy of surface hopping simulations in complex systems from quantum dynamics data.

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.

Detailed Balance and Independent Electron Surface-Hopping Method: The Importance of Decoherence and Correct Calculation of Diabatic Populations

We benchmark and improve the independent electron surface-hopping (IESH) method developed by J. C. Tully's group for nonadiabatic simulations near metal surfaces. We have incorporated decoherence within the IESH method as well as implemented a scheme for the accurate calculation of diabatic populations. We benchmark the original IESH method with the above inclusions for a model system to calculate rate constants and long-time populations. The original IESH method fails to capture the detailed balance for some of the parameters, which is corrected with the inclusion of decoherence and accurate calculation of diabatic populations. Total rate constants are well captured both within the original IESH method as well as within our modified IESH.

Independent Trajectory Mixed Quantum-Classical Approaches Based on the Exact Factorization

Mixed quantum-classical dynamics based on the exact factorization exploits the "derived" electron-nuclear correlation (ENC) term aiming for the description of quantum coherences. The ENC contains interactions between the phase of electronic states and nuclear quantum momenta which depend on the spatial shape of the nuclear density.The original surface hopping based on the exact factorization (SHXF) [\textit{J. Phys. Chem. Lett.} \textbf{2018}, \textit{9}, 1097] exploits frozen Gaussian functions to construct the nuclear density in the ENC term while the phase of electronic states is approximated as a fictitious nuclear momentum change.However, in reality, the width of nuclear wave packets varies in time depending on the shape of potential energy surfaces.In this work, we present a modified SHXF approach and a newly-developed Ehrenfest dynamics based on the exact factorization (EhXF) with time-dependent Gaussian functions and phases by enforcing total energy conservation.We perform numerical tests for various one-dimensional two-state model Hamiltonians.Overall, the time-dependent width of Gaussian functions and the energy conserving phase show a reliable decoherence compared to the original frozen Gaussian-based SHXF and the exact quantum mechanical calculation.Especially, the energy conserving phase is crucial for EhXF to reproduce the correct quantum dynamics.

Multilayer Subsystem Surface Hopping Method for Large-Scale Nonadiabatic Dynamics Simulation with Hundreds of Thousands of States

We present a multilayer subsystem surface hopping (MSSH) method to deal with nonadiabatic dynamics in large-scale systems. A small subsystem instead of the full system is adopted for surface hopping and is updated on-the-fly to achieve a reliable description of important adiabatic states and the wave function evolution. Additional subsystems for molecular dynamics and statistical description are introduced to further improve the simulation reliability. The global flux hopping probabilities with optimal state assignments are utilized to treat the complex surface crossings. As demonstrated in a series of one- and two-dimensional Holstein models with up to hundreds of thousands of states, MSSH shows weak parameter dependence in all investigated systems. Especially, the computational costs are reduced by 2-6 orders of magnitude compared to traditional surface hopping simulations in full systems, and size-independent results are achieved with a large time-step size of 2-5 fs. The new method is compatible with different decoherence correction strategies and achieves a much better balance between efficiency and reliability, thus promising for applications in general charge and exciton dynamics simulations.

A Mixed Deterministic-Stochastic Algorithm of the Branching Corrected Mean Field Method for Nonadiabatic Dynamics

We present a new algorithm of the branching corrected mean field (BCMF) method for nonadiabatic dynamics [J. Xu and L. Wang, J. Phys. Chem. Lett. 11, 8283 (2020)], which combines the key advantages of the two existed algorithms, i.e., the deterministic BCMF algorithm based on weights of trajectory branches (BCMF-w) and the stochastic BCMF algorithm with random collapse of the electronic wavefunction (BCMF-s). The resulting mixed deterministic-stochastic BCMF algorithm (BCMF-ws) is benchmarked in a series of standard scattering problems with potential wells on the excited-state surfaces, which are common in realistic systems. In all investigated cases, BCMF-ws holds the same high accuracy while the computational time is reduced about two orders of magnitude compared to the original BCMF-w and BCMF-s algorithms, thus promising for nonadiabatic dynamics simulations of general systems.

All-Atom Nonadiabatic Dynamics Simulation of Hybrid Graphene Nanoribbons Based on Wannier Analysis and Machine Learning

Trajectory surface hopping combined with ab initio electronic structure calculations is a popular and powerful approach for on-the-fly nonadiabatic dynamics simulations. For large systems, however, this remains a significant challenge because of the unaffordable computational cost of large-scale electronic structure calculations. Here, we present an efficient divide-and-conquer approach to construct the system Hamiltonian based on Wannier analysis and machine learning. In detail, the large system under investigation is first decomposed into small building blocks, and then all possible segments formed by building blocks within a cutoff distance are found out. Ab initio molecular dynamics is carried out to generate a sequence of geometries for each equivalent segment with periodicity. The Hamiltonian matrices in the maximum localized Wannier function (MLWF) basis are obtained for all geometries and utilized to train artificial neural networks (ANNs) for the structure-dependent Hamiltonian elements. Taking advantage of the orthogonality and spatial locality of MLWFs, the one-electron Hamiltonian of a large system at arbitrary geometry can be directly constructed by the trained ANNs. As demonstrations, we study charge transport in a zigzag graphene nanoribbon (GNR), a coved GNR, and a series of hybrid GNRs with a state-of-the-art surface hopping method. The interplay between delocalized and localized states is found to determine the electron dynamics in hybrid GNRs. Our approach has successfully studied GNRs with >10 000 atoms, paving the way for efficient and reliable all-atom nonadiabatic dynamics simulation of general systems.

Branching and phase corrected surface hopping: A benchmark of nonadiabatic dynamics in multilevel systems

Since the seminal work of Tully [J. Chem. Phys. 93, 1061 (1990)], two-level scattering models have been extensively adopted as the standard benchmark systems to assess the performance of different trajectory surface hopping methods for nonadiabatic dynamics simulations. Here, we extend the branching and phase corrections to multilevel systems and combine them with both the traditional fewest switches surface hopping (FSSH) and its variant global flux surface hopping (GFSH) algorithms. To get a comprehensive evaluation of the proposed methods, we construct a series of more challenging and diverse three-level and four-level scattering models and use exact quantum solutions as references. Encouragingly, both FSSH and GFSH with the branching and phase corrections produce excellent and nearly identical results in all investigated systems, indicating that the new surface hopping methods are robust to describe multilevel problems and the reliability is insensitive to the definition of self-consistent hopping probabilities in the adiabatic representation. Furthermore, the branching correction is found to be especially important when dealing with strongly repulsive potential energy surfaces, which are common in realistic systems, thus promising for general applications.

Nonadiabatic dynamics simulation of photoinduced ring-opening reaction of 2(5H)-thiophenone with internal conversion and intersystem crossing

In the present work, the quantum trajectory mean-field approach, which is able to overcome the overcoherence problem, was generalized to simulate internal conversion and intersystem crossing processes simultaneously. The photoinduced ring-opening and subsequent rearrangement reactions of isolated 2(5H)-thiophenone were studied based on geometry optimizations on critical structures and nonadiabatic dynamics simulations using this method. Upon 267 nm irradiation, the molecule is initially populated in the 1ππ* state. After a sudden rupture of one C-S bond within 100 fs in this state, the lowest two singlet excited states and the lowest two triplet excited states become quasi-degenerated, and then the intersystem crossing processes between singlet and triplet states accompanied by rearrangement reactions can be observed several times. Compared with our previous nonadiabatic simulations in the absence of intersystem crossing (ChemPhotoChem, 2019, 3, 897-906), some new nonadiabatic relaxation pathways involving triplet states and different ring-opening products were identified. The present work provides new mechanistic insights into the photoinduced ring-opening of thio-substituted heterocyclic molecules and reveals the importance of nonadiabatic dynamics simulation that is able to deal with multiple electronic states with different spin multiplicities.

Evaluation of mixed quantum-classical molecular dynamics on cis-azobenzene photoisomerization

The quantitative prediction of nonadiabatic transitions between different electronic states is important to understand ultrafast processes in photochemistry. A variety of mixed quantum-classical molecular dynamics methods such as surface hopping and Ehrenfest mean-field have been developed. However, how to choose an appropriate one from a wide diversity of dynamics algorithms to study a realistic photochemical process is still unclear. In this work, we implemented 30 combinations of different mixed quantum-classical dynamics methods, including 24 surface hopping models with 8 decoherence corrections and 3 momentum rescaling strategies as well as 6 mean-field models. Then we performed numerical investigations by simulating the photoisomerization of cis-azobenzene combined with on-the-fly electronic structure calculations. Predictions of the S1 lifetime and the quantum yield of the photoproduct using different models are distinct. Surface hopping is more robust than mean-field in our test system. Moreover, the choice of momentum rescaling methods in surface hopping brings more significant changes than decoherence corrections, while a large discrepancy between simulation results with different mean-field algorithms has been observed.

Branching Corrected Mean Field Method for Nonadiabatic Dynamics

When nonadiabatic dynamics are described on the basis of trajectories, severe trajectory branching occurs when the nuclear wave packets on some potential energy surfaces are reflected while those on the remaining surfaces are not. As a result, the traditional Ehrenfest mean field (EMF) approximation breaks down. In this study, two versions of the branching corrected mean field (BCMF) method are proposed. Namely, when trajectory branching is identified, BCMF stochastically selects either the reflected or the nonreflected group to build the new mean field trajectory or splits the mean field trajectory into two new trajectories with the corresponding weights. As benchmarked in six standard model systems and an extensive model base with two hundred diverse scattering models, BCMF significantly improves the accuracy while retaining the high efficiency of the traditional EMF. In fact, BCMF closely reproduces the exact quantum dynamics in all investigated systems, thus highlighting the essential role of branching correction in nonadiabatic dynamics simulations of general systems.

Mixed quantum-classical approach to model non-adiabatic electron-nuclear dynamics: Detailed balance and improved surface hopping method

We develop a density matrix formalism to describe coupled electron-nuclear dynamics. To this end, we introduce an effective Hamiltonian formalism that describes electronic transitions and small (quantum) nuclear fluctuations along a classical trajectory of the nuclei. Using this Hamiltonian, we derive equations of motion for the electronic occupation numbers and for the nuclear coordinates and momenta. We show that, in the limit, when the number of nuclear degrees of freedom coupled to a given electronic transition is sufficiently high (i.e., the strong decoherence limit), the equations of motion for the electronic occupation numbers become Markovian. Furthermore, the transition rates in these (rate) equations are asymmetric with respect to the lower-to-higher energy transitions and vice versa. In thermal equilibrium, such asymmetry corresponds to the detailed balance condition. We also study the equations for the electronic occupations in the non-Markovian regime and develop a surface hopping algorithm based on our formalism. To treat the decoherence effects, we introduce additional "virtual" nuclear wave packets whose interference with the "real" (physical) wave packets leads to the reduction in coupling between the electronic states (i.e., decoherence) as well as to the phase shifts that improve the accuracy of the numerical approach. Remarkably, the same phase shifts lead to the detailed balance condition in the strong decoherence limit.