SHARC meets TEQUILA: mixed quantum-classical dynamics on a quantum computer using a hybrid quantum-classical algorithm

Abstract

Recent developments in quantum computing are highly promising, particularly in the realm of quantum chemistry. Due to the noisy nature of currently available quantum hardware, hybrid quantum-classical algorithms have emerged as a reliable option for near-term simulations. Mixed quantum-classical dynamics methods effectively capture nonadiabatic effects by integrating classical nuclear dynamics with quantum chemical computations of the electronic properties. However, these methods face challenges due to the high computational cost of the quantum chemistry part. To mitigate the computational demand, we propose a method where the required electronic properties are computed through a hybrid quantum-classical approach that combines classical and quantum hardware. This framework employs the variational quantum eigensolver and variational quantum deflation algorithms to obtain ground and excited state energies, gradients, nonadiabatic coupling vectors, and transition dipole moments. These quantities are used to propagate the nonadiabatic molecular dynamics using the Tully's fewest switches surface hopping method, although the implementation is also compatible with other molecular dynamics approaches. The approach, implemented by integrating the molecular dynamics program package SHARC with the TEQUILA quantum computing framework, is validated by studying the cis–trans photoisomerization of methanimine and the electronic relaxation of ethylene. The results show qualitatively accurate molecular dynamics that align with experimental findings and other computational studies. This work is expected to mark a significant step towards achieving a “quantum advantage” for realistic chemical simulations.

---

1 Introduction

As suggested by Feynman¹ in the 1980s, simulating quantum mechanical systems on a quantum computer should be more efficient than running the same computation on a classical computer. Especially quantum chemical simulations have emerged as particular promising applications.^2–5 Current quantum computing hardware classify as “noisy intermediate scale quantum technology” (NISQ) devices.^6,7 The limited amount of available qubits and the still open problem of efficient error correction leaves the complete simulation of large and complex physical and chemical systems a distant goal. Hybrid quantum-classical algorithms, where only parts of the problem are solved on quantum devices⁸ are expected to be crucial in near-term quantum simulations,⁶

Over the past decade the variational quantum eigensolver (VQE) has become the de-facto standard hybrid quantum-classical algorithm for quantum chemical calculations.^9–11 Despite its initial design for computing electronic ground states, different variations for calculating electronically excited states have been developed. Especially penalty-based optimization^12–15 and quantum equations of motion (qEOM)¹⁶ approaches have shown promising results.

Mixed quantum-classical dynamics methods, which simulate nonadiabatic processes by combining classical trajectory propagation for the nuclear dynamics with a quantum mechanical treatment of the electrons, are effective in capturing nonadiabatic effects through a feedback algorithm that integrates the electronic and nuclear subsystems. Using a classical treatment of the nuclei allows for a favorable scaling with the system size. The significant bottleneck in these simulations is the computation of the electronic properties. The most common approach in mixed quantum-classical dynamics is to compute the electronic properties on-the-fly, i.e., at each timestep of the dynamics. Advantageously, this does not require pre-computed multidimensional potential energy surfaces. Still the computational cost of these dynamics simulations is heavily dependent on the accuracy of the electronic property calculations, making it necessary to find a compromise between accuracy and feasibility.

Moving portions of the quantum chemical computations for the electronic structure onto quantum computers should accelerate the overall computational effort. In this work, we utilize hybrid quantum-classical algorithms to compute electronic properties on-the-fly during mixed quantum-classical dynamics. Recently, some studies^17,18 have demonstrated the ability of hybrid quantum-classical algorithms to accurately predict conical intersections, which is relevant to perform nonadiabatic dynamics.

While previous studies have explored the application of hybrid quantum-classical algorithms in the context of excited state dynamics, they are often limited to toy models and very small molecular systems.^19–24 In contrast, we focus on a hybrid quantum-classical framework able to perform nonadiabatic dynamics for polyatomic molecules with complex and interesting photochemistry.

Although our approach is implemented using Tully's fewest switches surface-hopping (SH) method,²⁵ it can be easily extended to other molecular dynamical approaches.

The ground state electronic properties are computed using a standard VQE algorithm, while excited state properties are obtained through penalty-based optimization, which ensure the orthogonality of the wavefunction with previously calculated electronic states, thus allowing to subsequently compute multiple electronically excited states. Our approach is implemented by interfacing our molecular dynamics program package SHARC^26,27 with the TEQUILA quantum computing framework.²⁸

Recently, Tavernelli and co-workers employed the quantum subspace expansion and quantum equation-of-motion algorithms to simulate the SH dynamics of a hydrogen atom colliding with a hydrogen molecule, a system comprising three nuclei and three electrons.²² In the subspace expansion framework, excited states are represented as linear combinations of intermediary states, typically derived from the approximate ground-state wave function. While effective, the use of truncated subspaces is not strictly variational for excited states, which can sometimes limit accuracy. In contrast, we have employed a penalty-based method (see Section 2.2). This is not only simpler to implement, but also eliminates the dependence on the number of excited states required for accuracy. In our approach, each additional state is calculated variationally through a penalty function, enforcing orthogonality to previously computed states and maintaining precision without the limitations of truncated subspaces.

As we are bridging both the nonadiabatic dynamics and quantum information communities, the next section explains how electronic properties are computed using the VQE algorithm. In the following, we provide a detailed outline of the dynamics involved and illustrate how these electronic properties, obtained through the VQE, enter the nuclear dynamical procedure. Finally, we showcase our approach by investigating the nonadiabatic dynamics of two systems: the cis–trans photoisomerization of methanimine and the ultrafast electronic relaxation of ethylene.

2 Methods

2.1 Ground state energies and gradients

We utilized the hybrid quantum-classical algorithm, the VQE, to compute the ground state energies and gradients during the dynamics. The VQE was first introduced by Peruzzo et al.⁹ and subsequently extended later in ref. 10. It is based on the variational principle, which involves optimizing an upper bound for the lowest possible expectation value of an observable with respect to the trial wavefunction.¹¹ Specifically, given a Hamiltonian and a trial wavefunction Ψ(θ), the ground state energy associated with this Hamiltonian, E_gs, is bounded by:where the trial wavefunction Ψ(θ) is assumed to be normalized, and θ is the set of variational parameters. Therefore, the objective of the VQE is to find the set of θ that minimizes the expectation value of . The trial wavefunction is constructed using a quantum circuit, known as the ansatz. This involves applying a series of quantum gates (unitary operations) to a reference state Ψ₀. Thus, the trial wavefunction Ψ(θ) can be represented as U(θ)|Ψ₀〉, where U(θ) denotes the ansatz.

The VQE algorithm has shown to be an interesting approach for solving quantum chemistry problems. In this context, the molecular Hamiltonian in second quantization can be written as:

where h_rs and g_pqrs denote the one- and two-electron integrals, respectively, and are commonly obtained from a Hartree–Fock (HF) calculation, which is usually performed on a classical device. The collective vector of nuclear coordinates R = (R₁, R₂, …, R_N) of N nuclei in R^3N simply parametrizes the electronic Hamiltonian. The position vector of a single nucleus I ∈ {1, …, N} is denoted by R_I = (R_Ix, R_Iy, R_Iz). The operators and â_r represent the fermionic creation and annihilation operators for electrons in HF spin–orbitals. The indices r, s, t, u are used to label general (occupied or virtual) molecular orbitals. The term E_NN(R) describes the nuclear repulsion energy.

Observables suitable for direct measurements on a quantum device are tensor products of spin operators (Pauli operators). To achieve this, the Hamiltonian of eqn (2) can be transformed into a qubit operator by applying, for example, the Jordan–Wigner transformation,²⁹ so that the molecular Hamiltonian can be expressed as:

where w_α is the set of weights and the number of Pauli strings P̂_α in the Hamiltonian.

Therefore, the expectation value of the molecular Hamiltonian with the VQE can be calculated with:

where the hybrid nature of the VQE algorithm becomes clearly apparent: each term E_Pα = 〈Ψ₀|U^†(θ)P̂_αU(θ)|Ψ₀〉 corresponds to the expectation value of a Pauli string P̂_α and can be computed on a quantum device, while the summation and minimization is computed on a conventional computer. The reference wavefunction |Ψ₀〉 is often the HF wavefunction. The evaluation of the expectation value in eqn (3) often becomes the computational bottleneck of the VQE procedure, as the number of Pauli strings in a naive decomposition scales quartic with respect to the number of orbitals. There exist however methods to group the Pauli strings into commuting cliques^30–32 that mitigates this issue – the number of diagonal bases in the ethylene molecule used below is for example reduced from 34 (native Pauli bases) to just 3 (commuting Pauli cliques) using the method of ref. 31. For a full treatment (6-31G basis), the reduction would be from 73 089 to 552.

The quality of the energy obtained from eqn (4) also depends on the ansatz. Several ansätze have been proposed in recent years for VQE applications in quantum chemistry. Among the most commonly used is the Unitary Coupled Cluster (UCC) approach and its extensions.^33–36 Other popular choices include the Hardware Efficient Ansatz (HEA),³⁷ and the anti-Hermitian contracted Schrödinger equation (ACSE) ansatz, which has been effectively used in contracted quantum eigensolvers (CQE) on quantum computers.^38–40 In this work, we have adopted a special type of UCC, termed k-UpCCGSD.⁴¹ This ansatz presents a good compromise between accuracy and cost, showing better scaling than UCCSD and UGCCSD.

The derivatives of the energy with respect to the nuclear coordinates for a given nucleus I, (energy gradients), where ξ ∈ {x, y, z}, are needed to calculate the forces acting on the nuclei. The total derivative is given explicitly by:

with . The first term in eqn (5) is the Hellmann–Feynman force, while the second and third terms are the Pulay forces. In our implementation, we approximated the energy gradients by considering only the Hellmann–Feynman term:where the gradients of the molecular Hamiltonian were calculated using the central difference formula with a shift parameter of 0.001 angstroms.

We note that this approximation has also been used in other studies that disregard Pulay forces.^19–22 Already in 2010, Foley and Mazziotti⁴² demonstrated efficient geometry optimization and frequency analysis by using only the Hellmann–Feynman component for the ACSE. Pulay forces are considered essential for correcting errors arising from basis set incompleteness.⁴³ Consequently, a complete basis set would result in net zero Pulay forces. Since we did not use a complete basis set, we assessed the accuracy of neglecting Pulay forces by comparing gradients calculated solely using the Hellmann–Feynman term with those obtained using the central difference formula of the energy (which includes both Pulay and Hellmann–Feynman terms) across various geometries. Generally, the gradients from both methods were nearly identical, indicating that neglecting Pulay forces is valid in this context. Further details and discussions on these results are provided in Section S1 of the ESI.†

It is also worth mentioning that another alternative, which allows forces to be computed solely using the Hellmann–Feynman term, is to use a perturbation-dependent basis set.⁴⁴

2.2 Excited state energies and gradients

The excited states can be calculated similarly to the method used by the VQE, but instead of minimizing the energy, a cost function is minimized.^12–15 This method is known as variational quantum deflation (VQD). The definition of the cost function depends on the target state and the target system. For instance, the energy of the first excited state can be calculated using the cost function:Here, Ψ_gs is the ground state wavefunction that was previously computed with VQE, and λ is a hyperparameter that must be specified before the VQD calculation. The second term in eqn (7) enforces the constraint of searching the subspace orthogonal to the ground state. Following the ref. 45, we can express the overlap punishment as another expectation value

|〈Ψ(θ)|Ψ_gs〉|² = 〈Φ(θ)|P₀|Φ(θ)〉 (8)

with |Φ(θ)〉 = U(θ)^†|Ψ(θ)〉 and P₀ = |0〉〈0| with the all-qubit zero state |0〉. The parameter λ should be sufficiently large, ideally larger than the energy difference between the ground and the excited state. However, if λ is too large, it could lead to the selection of an undesired higher excited state.⁴⁶ A good initial guess for λ could be the magnitude of the ground state energy obtained with VQE – which is sufficient if the excited state is bound (i.e. has negative energy)⁴⁵ – and this is the λ used in this work. In addition, the ansatz adopted in this work, unlike the HEA, prevents spin contamination since the number of electrons and the spin state are conserved. Therefore, no additional considerations, such as including penalty terms in the cost function (eqn (7)) are needed.

We need to stress that the VQD faces significant challenges due to its sequential computation requirement. Each eigenstate computation depends on the accurate deflation of previously found states, leading to increased computational overhead and potential optimization difficulties. This sequential nature can complicate the optimization process, making it harder to accurately identify and optimize higher excited states. Consequently, VQD may miss some states or converge to suboptimal solutions, particularly in more complex systems or under noisy quantum conditions.

The gradients of the excited state energies are obtained similarly to those for the ground state (see eqn (6)). Once we have the gradient of the molecular Hamiltonian, , the gradients of all the states can be easily calculated as the expectation value of this derivative for the wavefunctions obtained from the VQE or VQD algorithms. Therefore, the computational cost for the computation of the gradients does not increase significantly with the number of states.

2.3 Surface hopping dynamics

The time evolution of a (non-relativistic) molecular system is determined by the time-dependent Schrödinger equation (TDSE), . Unfortunately, solving the TDSE and simulating quantum dynamics in molecules remains a fundamentally difficult problem due to the non-local nature of quantum mechanics and the associated exponential scaling of the computational effort with the number of degrees of freedom. One popular and effective approximation is to describe the nuclear motion using classical mechanics instead of quantum mechanics. This approach is motivated by the favorable scaling of classical mechanics with system size. The classical approximation yields reasonable results in cases where quantum mechanical effects, such as tunneling or interference, are negligible, and the energetic spacing between quantum levels is small compared to the kinetic energy. Since nuclei are much heavier than electrons, the spacing between quantum levels is usually much smaller for nuclear degrees of freedom than for electronic ones. Therefore, it is plausible to treat the nuclei with classical mechanics while treating the electrons with quantum mechanics. Based on this assumption, numerous quantum-classical methods have been developed over the past decades to address nonadiabatic processes. Among these methods, the SH scheme has emerged as one of the most widely used approaches.^47–50 In this work, we used the Tully's fewest switches SH,²⁵ but the approach presented here can be easily extended to others mixed quantum-classical dynamics approaches. The basic assumption is that during nonadiabatic dynamics, the nuclei move adiabatically for most of the time and only undergo nonadiabatic transitions for relatively short periods and in relatively small regions of the configuration space. Hence, it is pragmatically proposed that one could approximate the nonadiabatic transitions by instantaneous “hops” between adiabatic potential energy surfaces (PESs). In the SH approach, many trajectories are simulated independently, generating a statistical ensemble. The fraction of trajectories for each PES is used to mimic the population of each quantum state in realistic dynamical processes. In SH dynamics, the time evolution of the nuclei, i.e., the classical degrees of freedom R, is performed by integrating Newton's equation with the potential being in a given adiabatic surface E_m(R). Here, Newton's equation is solved using the velocity-Verlet algorithm.⁵¹ In this algorithm, the nuclear positions R_I,ξ and velocities Ṙ_I,ξ are updated from a time step, t, to the next step, t + Δt, with the following equations:where M_I is the atomic mass of a given nucleus I, and the forces acting on the nuclei, F_I,ξ, are the negative derivative of the electronic energy with respect to the nuclear coordinates, .

Besides the nuclei, the electronic wavefunction must be propagated as well. The electronic wavefunction is expressed as a sum over basis states:

where m indexes the basis states |ψ_m〉 (eigenstates of the electronic Hamiltonian), a_m(t) are time-dependent coefficients, and γ_m(t) represents the phase: with E_m(R) being the adiabatic energy evaluated at the nuclear configuration R.

Incorporating Ψ_el(t) into the TDSE yields the equation for the time evolution of the coefficients a_m(t):

where σ_ml represents the time-derivative nonadiabatic coupling:with v being the nuclear velocity, and g_ml the first-order nonadiabatic coupling (NAC) vector, defined as:Within the VQE/VQD framework, the NAC vectors can be approximated aswhere Ψ(θ_m) and Ψ(θ_l) represent the wavefunctions optimized by VQE or VQD for states m and l, respectively, with corresponding energies E(θ_m) and E(θ_l). Once the energies and gradients are obtained, evaluating the NAC vectors is straightforward and present a similar computational cost to that of computing the gradients (see eqn (6)). A similar methodology for computing first-order NAC vectors was proposed in ref. 24.

However, NAC vectors can change rapidly, often necessitating very small time steps to avoid numerical issues, especially near narrow peaks found in NACs near strong interaction regions such as avoided crossings and conical intersections. Additionally, crossings between uncoupled states (“trivial crossings”) are invariably missed regardless of time step size.⁵² To circumvent these challenges, we employed the local diabatization scheme for integrating electronic coefficients.^53,54 In this scheme, the NAC vectors are never explicitly computed. Instead, the wavefunction overlap matrix at two different time steps is required. The basic idea of the local diabatization is to resort to a “locally diabatic” representation, i.e., to a set of electronic states that are specifically diabatic along the nuclear trajectory under consideration. By definition, the diabatic basis also spans the internal subspace and is connected with the adiabatic one by a unitary transformation:

|η〉 = |ψ〉T (14)

where |η〉 and |ψ〉 are the wavefunctions in diabatic and adiabatic basis, respectively.

If H is the Hamiltonian in the diabatic basis (H_kl = 〈η_k|Ĥ_el|η_l〉) and E the diagonal matrix of the electronic energies, we have

HT = TE. (15)

The diabatic expansion of the time-dependent wavefunction is

with and

d = Ta (17)

where d and a are the matrices of the diabatic and adiabatic coefficients, respectively. Here, H_mm are the diagonal elements of the Ĥ_el matrix in the diabatic basis, i.e., H_mm = 〈η_m|Ĥ_el|η_m〉.

The diabatic basis in the local diabatization scheme is redefined at each time step. At the beginning of the time step, the transformation matrix is chosen to be the identity matrix, thus η ≡ ψ. By choosing the diabatic states to be constant within the time step considered, the dynamic couplings 〈η_k|∂/∂t|η_l〉 are eliminated. However, the coupling vanishes only along the advancement coordinate identified by the velocity vector Ṙ: in this sense the η are “locally” diabatic” for the given trajectory. Then, by inserting the electronic wavefunction expansion (16) into the TDSE, the time evolution of the diabatic coefficients is obtained:

Within the local diabatization algorithm, the transformation matrix T is obtained by Löwdin orthonormalization of the wavefunction overlap matrix

S_ml (t + Δt) = 〈ψ_m(t)|ψ_l(t + Δt)〉. (19)

Hence, from T, one obtains H(t + Δt). The matrix H at intermediate times in the interval [t, t + Δt] is attained by interpolation from H(t) to H(t + Δt) (note that H(t) = E(t)). The essential point is that, unlike the NACs, the diabatic quantities (such as H) depend smoothly on the nuclear coordinates and can be easily and effectively interpolated. Then, the diabatic coefficients are obtained from eqn (18). Finally, the adiabatic expansion coefficients a(t + Δt) are computed through the inverse of eqn (17).

In the present context, ψ_m represents the wavefunctions optimized by VQE or VQD for states m and l, denoted as Ψ(θ_m) and Ψ(θ_l), respectively. However, these wavefunctions correspond to different time steps, meaning that the molecular orbitals in the wavefunctions vary over time. However, to compute the overlaps in eqn (19), we must first convert the wavefunction into a “configuration interaction type” wavefunction. This process requires providing a list of determinants, their corresponding coefficients, and the molecular orbitals for both the bra and ket states. These overlaps are then computed using a classical algorithm.

The local diabatization scheme allows using larger time steps and naturally mitigates the trivial crossing issue.

2.3.1 Initial conditions preparation. We showed how in SH the nuclear and electronic properties (adiabatic energies, energy's gradients and NACs) can be propagated from a time step t to the next step, t + Δt. However, to start an SH trajectory, we need to provide: (i) the initial nuclear positions R; (ii) initial nuclear velocities v; (iii) the initial electronic state m; and (iv) the initial matrix of the electronic wavefunction coefficients a. This set of values is called the initial conditions.

The electronic wavefunction coefficients a and the initial state m depend on the excitation process. There are two techniques that are often used to sample R and v. The first one is called Wigner sampling and it computes an approximate phase space probability distribution function representing the ground vibrational state and then stochastically draws samples from this distribution. The advantage of Wigner sampling is that quantum effects like zero-point energy are adequately considered. However, it is often challenging to find an appropriate nuclear wavefunction for polyatomic molecules; thus, one usually resorts to the harmonic approximation, which works well for small and stiff molecules.^55,56 An alternative to Wigner sampling is to run a long dynamics simulation in the ground state and randomly pick R and v from this trajectory. This approach does not represent quantum effects like zero-point energy well but works effectively for large, polyatomic systems with anharmonic or nonlinear modes and multiple local minima in the ground state PES.^57–59 While the latter approach is easier to implement, as it only requires energies and gradients-which we have already demonstrated how to compute-we decided to use Wigner sampling, as we believe this method is more appropriate for the types of systems studied in this work. For Wigner sampling, one needs to provide the normal modes and their corresponding frequencies. This requires calculating the second derivatives of the energy with respect to the nuclear coordinates. Since obtaining the analytical solution for this derivative is not straightforward, and since we only need to compute these derivatives at a single point-specifically, the minimum of the ground state, we opted for a fully numerical approach using central finite differences to calculate second derivatives.

As we have sampled the quantities R and v, it is possible to find the electronic wavefunction coefficients a and the initial state m. In the simplest case, it is possible to set a_m(0) = δ_ml, with l being one of the excited states. This l can be simply defined by the user (e.g., all trajectories start in l = 3, i.e., the S₂ state).

A more appropriate and popular approach is to perform a single point calculation for each sampled R and compute a selection probability based on the obtained excitation energies and oscillator strengths of each state

where f_m is the oscillator strength of state m, and p_norm is chosen arbitrarily. With these probabilities, one can stochastically select initial states. Additionally, some restrictions on |E_m(R) − E_gs(R)| are imposed to consider only a small excitation energy window. Therefore, our implementation also includes the calculation of transition dipole moments within the VQE/VQD framework. These calculations are necessary to determine oscillator strengths and to select the initial conditions accordingly.

2.4 Dynamics workflow

In our implementation, the computation of the electronic properties such as energies, gradients, and transition dipole moments are carried out using the TEQUILA quantum computing software²⁸ using the Qulacs⁶⁰ quantum circuit simulation library. The overlap matrix between consecutive time steps (eqn (19)) is computed using the WFoverlap program,⁶¹ and the SH dynamics are executed using the SHARC package.^26,27 After selecting the initial conditions, the workflow for each SH trajectory in the VQE/VQD framework is as follows:

(1) Calculate energies according to eqn (4) and (7).

(2) Compute gradients using eqn (6).

(3) Determine the overlap matrix between consecutive time steps (eqn (19)).

(4) Provide energies, gradients, and the wavefunction overlap matrix to the driver responsible for nuclear dynamics—the SHARC software in this case. This driver integrates the electronic TDSE, evaluates hopping probabilities, and computes new nuclear coordinates and velocities.

(5) Repeat steps 1–4 until reaching the desired trajectory endpoint. Fig. 1 illustrates a schematic representation of our workflow, highlighting the components performed on classical and quantum devices.

Fig. 1 General workflow used to perform surface hopping dynamics in combination with the VQE/VQD algorithm. All parts of the simulations running on classical hardware (CPU) are depicted on the left hand side of the picture and shaded in blue. The preparation of the electronic wavefunction and the evaluation of the ground and excited state energies as well as the gradients and (optionally) transition dipole moments and NACs are carried out on the quantum device (QPU), shown on the right hand side of the picture, shaded in red. These data are then used by the classical hardware to propagate the nuclear dynamics (elements highlighted in white on the left panel).

3 Applications of surface-hopping dynamics on quantum computers

In order to validate our SH implementation on quantum computers, we applied it to two case studies: the cis–trans photoisomerization of methanimine (HN=CH₂) and the electronic relaxation of ethylene (CH₂=CH₂). Both molecular systems are shown in Fig. 2. The small molecular size, ultrafast dynamics through a conical intersection, and massive conformational flexibility make these molecules suitable for ensure that our approach works as expected. Moreover, a large number of nonadiabatic dynamical studies have been carried out on these systems, employing various dynamical methods combined with different levels of theory to describe the PESs,^62–80 which serve as reference. Both systems are also ideal to test our algorithm, as a proper description of the electronic structure problem only requires a small number of active orbitals. This is particularly important as the number of qubits required for the ansatz used is proportional to the number of active spin–orbitals.

Fig. 2 Left: Methanimine (HN=CH₂). Right: Ethylene (CH₂=CH₂). Carbon (C) atoms in gray, nitrogen (N) atom in blue and hydrogen (H) atoms in white.

3.1 Test case HN=CH₂

3.1.1 Computational details. The VQE/VQD electronic energies and wavefunctions of HN=CH₂ were obtained using an active space of 4 electrons and 3 orbitals (6 spin–orbitals), requiring 6 qubits to construct the ansatz. The active space consisted of the n, π, and π* orbitals, warrantying that both nπ* and ππ* excitations were considered. A HF/6-31G wavefunction was used as a reference, calculated with PySCF.⁸¹ The Jordan–Wigner mapping²⁹ was applied to generate the qubit Hamiltonian, and both VQE and VQD used the BFGS optimizer to minimize the parameters θ. In the first time step, the parameters θ were set to zero. However, for subsequent time steps, the initial θ values were those obtained from the optimization of the previous time step. This approach serves two purposes: first, it is expected that fewer interactions will be needed to optimize θ for the new geometry rendering the computations faster; second, and more importantly, it ensures that the VQD algorithm tracks the same set of states throughout the dynamics. For the first excited state, 6 angles were optimized, requiring an average over a single trajectory of 8.12 iterations (optimization steps) per nuclear time step (with a time step of 0.5 fs). For the second excited state, 13 angles were utilized, but only 7 were optimized, necessitating an average of 13.87 iterations. For the third excited state, 19 angles were involved, with 6 optimized, averaging 20.04 iterations per nuclear time step. Each nuclear time step involved a new calculation, resulting in varying numbers of interactions during the dynamics. As mentioned previously, the optimized angles from the previous time step were used as initial guesses for the next step, which reduced the number of interactions needed. Additionally, as expected, reducing the nuclear time step led to fewer iterations; for example, when the time step was reduced to 0.2 fs, the number of iterations was approximately halved.

The distribution of initial coordinates and velocities was generated from a Wigner distribution, from which the spectrum shown in Fig. 3 is obtained. The starting conditions (geometries, velocities, and initial state) were selected according to the excitation window ranging from 5.5 eV to 6.0 eV, which corresponds to the excitation to the S₁ state (nπ* excitation). The lowest three lying electronic states (S₀–S₂) were considered during the SH simulations. A total of 150 initial conditions were extracted from a nuclear ensemble of 5000 geometries. In total, 141 were used to analyze the dynamics; 9 trajectories were discarded because the total energy was not conserved by the end of the dynamics due to changing the orbitals in the active space, resulting in inconsistent PESs or numerical instabilities in the gradient computation. While the orbital rotation problem is typically resolved by increasing the active space size, here we keep the small active space consisting of only three orbitals, as only very few trajectories are affected by this problem.

Fig. 3 Simulated absorption spectrum of HN=CH₂ calculated using 5000 geometries. The inset shows the nπ* band. Contributions from each adiabatic state are indicated, with the dashed line representing the total spectrum.

The local diabatization algorithm was used for the integration of the electronic TDSE,^53,54 with an integration time step of 0.02 fs for the electronic degrees of freedom and 0.5 fs for the nuclear degrees of freedom. The Granucci–Persico energy-based decoherence correction with the standard 0.1 a.u. parameter was applied during the SH dynamics.⁸² Rescaling of the nuclear velocities after a hop was performed in the direction of the nuclear momentum. The SH trajectories were propagated up to 150 fs but, in order to save computational time, were stopped if reached the electronic ground state and stayed there for at least 20 fs.

3.1.2 Simulation of the excited-state dynamics. The simulated absorption spectrum, shown in Fig. 3, exhibits two peaks: a weaker one at lower energy corresponding to the forbidden nπ* excitations, and a stronger peak at higher energies corresponding to the ππ* excitation. The maximum of the nπ* state band (∼5.7 eV) qualitatively reproduces the vertical excitation obtained with previous results, e.g. with spin–flip time-dependent density functional theory (SF-TDDFT) (5.65 eV),⁷⁷ linear-response time-dependent density functional theory (LR-TDDFT) (4.84 eV and 4.92 eV),^76,79 configuration interaction (∼5 eV),⁸³ quantum Monte Carlo (5.32 eV),⁸⁴ and the experiment (∼4.96 eV).⁸⁵

HN=CH₂, the smallest unprotonated Schiff base, is a prototypical system for studying cis–trans photoisomerization around a double bond. The photophysical properties and associated nonadiabatic photoisomerization mechanism for this system have been rigorously studied using nonadiabatic ab initio molecular dynamics with LR-TDDFT,⁷⁶ SH dynamics with semiempirical potentials and SF-TDDFT,^77,78 and high-level “static” multireference configuration interaction calculations.⁸³ In the electronic ground state, S₀, HN=CH₂ adopts the planar geometry depicted in Fig. 2. The system contains five atoms and therefore nine internal degrees of freedom that are all free to relax during the excited state dynamics simulations. After vertical excitation to the lowest excited singlet state, S₁, it is known that the system rapidly relaxes toward the local energy minimum on the excited state PES and the N–H bond vector twists out of the molecular plane.^76–78 On approaching the orthogonal twist geometry, the system enters the conical intersection region, resulting in strong nonadiabatic coupling between the S₁ and S₀ states. After relaxation to the S₀ state, HN=CH₂ returns to a planar geometry, leading either to the photoisomerized product or to the regeneration of the reactant state.

Our simulations reproduce a similar behavior, as shown by the example trajectory in Fig. 4. Initially, in structure (a), the molecule presents a planar configuration. Structure (b) represents the hopping structure, where the N–H bond vector twists out of the molecular plane, followed by photoisomerization, leading to structure (c) in one of the last time steps. Other trajectories are shown in Section S3 (ESI).†

Fig. 4 Potential energy surfaces calculated with VQE/VQD for a typical reactive cis–trans trajectory of HN=CH₂, illustrating (a) the initial planar configuration, (b) the hopping structure with the N–H bond vector twisted out of the plane, and (c) the photoisomerized structure. The thin black line below 8 eV represents the total energy (kinetic + potential) over time and the circles represent the active adiabatic state.

Previous calculations using classical algorithms reported that the crossing point is characterized by an HNC angle and an HNCH dihedral of approximately 100°. For instance, Tavernelli et al. reported an HNC angle and an HNCH dihedral of ∼100°,⁷⁶ while Bonačić-Koutecký and Michl reported an HNC angle of ∼106.5° at the crossing point.⁸⁶ We calculated the average of these internal coordinates for all the hopping geometries between the S₀ and S₁ states. Our findings show an average HNC angle of 108.4° and an HNCH dihedral (dihedral among the atoms 1–2–3–4; see Fig. 2) of 100.5°, which is in very good agreement with previous studies carried out on purely classical computers. In Fig. 5 panel (b), we show the HNC angles and the HNCH dihedral at the hopping geometries between S₀–S₁ and S₁–S₂. We note that most of the hops take place when both angles and dihedrals are in the range of 80–120°. This confirms that our approach can accurately reproduce the conical intersection geometries and is stable even at strong nonadiabatic coupling regions.

Fig. 5 Panel (a) top: convolution of the HNC angle over time from 141 trajectories. Panel (a) bottom: convolution of the HNCH dihedral (dihedral among the atoms 1–2–3–4; see Fig. 2) over time from all trajectories. Panel (b): HNCH dihedral (dihedral among the atoms 1–2–3–4; see Fig. 2), in red, and HCN angle, in blue, at the S₁–S₀/S₀–S₁ hopping moment.

Fig. 5 panel (a) depicts the convolution of the HNC angle and the HNCH dihedral over time from all the trajectories considered. The HNC oscillates from 100–180° with a periodicity of 20 fs, while the HNCH dihedral oscillates from 0–180° with a periodicity of 40 fs. This periodicity diminishes as the trajectories start to decay to the ground state. Note that the strong coupling regime is achieved mainly when both internal coordinates are close to 100°.

In Fig. 6 the population dynamics of the adiabatic states averaged over all trajectories following nπ* excitation is shown. Initially, all populations are in the S₁ state, which has nπ* character, while the S₂ (ππ*) state remains unpopulated throughout the dynamics. This observation aligns with the energy gap between these states (see Fig. 3). From the beginning of the dynamics, until approximately 30–40 fs, only a few trajectories decay to the ground state S₀. Subsequently, an exponential decay to the ground state is observed.

Fig. 6 Time-resolved adiabatic state populations of HN=CH₂, averaged over 141 trajectories. The gray line represents the fit of the S₁ population according to the kinetic model described by eqn (21).

To quantify this decay, we fit the S₁ population P_S1 using a simple exponential decay model with a delay time t₀. This model is described by the equation:

where τ₁ is the decay constant, and the sum τ₁ + t₀ gives the lifetime of the S₁ population. Here, we used a delay time t₀ = 30.00 fs, which resulted in a decay constant τ₁ = 136 fs. Therefore, the lifetime of the S₁ population is 167 fs.

The observed population decay behavior is consistent with previous findings, e.g. those reported in ref. 78, where they used SH dynamics and semiempirical PESs, finding a similar “plateau” in the S₁ population until 25–30 fs, followed by an exponential decay to the ground state. Dynamical simulations using SF-TDDFT⁷⁷ showed a faster decay (∼58 fs) to the ground state compared to our observations. The differences can be attributed to minor variations in the PESs, that in return can have a substantial impact on the lifetimes. Additionally, the higher energy position of the nπ* band reported by us compared to the experimental spectrum⁸⁵ and previous TDDFT calculations⁷⁷ may contribute to the observed delay in the decay to the ground state.

Here, calculating the photoisomerization quantum yield is not straightforward because it can only be determined for trajectories that end in the ground state. However, we terminated the trajectories after 20 fs in the ground state, which is insufficient for the system to fully relax to the final product. As shown in Fig. 6, 62% of the population (corresponding to 88 trajectories) relaxed to the ground state by at the end of the propagation. Among these 88 trajectories, 30 exhibited an HNCH dihedral angle of ≤50°, suggesting that the final product is likely the non-photoisomerized form (we analysed the dihedral among the atoms 1–2–3–4; see Fig. 2), while 38 trajectories exhibited an HNCH dihedral angle ≥125°, indicating that the product most likely remains in the photoisomerized form. The remaining 20 trajectories, which also ended in the ground state, had dihedral angles between 50–125°, so we did not assign a final product to them. Considering only the 68 trajectories with definitive outcomes (30 cis and 38 trans), we estimate a “partial” photoisomerization quantum yield, Φ, of approximately 56% with a standard deviation of ±6%. This standard deviation is calculated according to the binomial standard deviation formula: . This result is consistent with the results obtained for SH dynamics using semiempirical PES (54%).⁷⁸ However, it is lower than the results from nonadiabatic Car–Parrinello molecular dynamics, which reported a value of 69%.⁸⁰

Finally, we investigated whether the VQD algorithm can always find states that are orthogonal to each other. Therefore, we calculated the overlap between the electronic states 〈Ψ(θ_m)|Ψ(θ_l)〉 over time. Note that this overlap is not the same as in eqn (19). That overlap pertains to states at different times, meaning they represent different sets of molecular orbitals and its computation is more sophisticated. In Section S2 (ESI†), we shows for four different trajectories, how the overlaps between states S₀–S₁, S₀–S₂, and S₁–S₂ evolve over time. We observe that the overlaps are generally very small, indicating that most of the time, the VQD converges to states that are orthogonal to each other. In a few instances during the dynamics, the overlap can reach values of the order of 10⁻³. However, these overlaps always involve the S₂ state, which does not participate in the dynamics. The overlaps between S₀–S₁, which are the states participating in the dynamics, remain very small, generally bellow to 1 × 10⁻⁵. From this, we conclude that the VQD is robust enough to reproduce the correct dynamical behaviour during the propagation.

3.2 Test case CH₂=CH₂

3.2.1 Computational details. A similar approach is used to perform the SH dynamics with slight differences.

The VQE/VQD electronic energies and wavefunctions were also obtained using an active space of 4 electrons and 3 orbitals (6 spin–orbitals and 6 qubits), which involves one π, one σ, and one π*. A HF/6-31G wavefunction was used as a reference, calculated with PySCF. The Jordan–Wigner mapping was applied to generate the qubit Hamiltonian, and both VQE and VQD used the BFGS optimizer to minimize the parameters θ. In the first time step, the parameters θ used were the optimal ones for HN=CH₂. In attempting to set them to zero, as we did at the first time step of HN=CH₂, the VQD did not find the state corresponding to the ππ* transitions (this state can be the S₁ or the S₂ state depending on the geometry). However, for subsequent time steps, the initial θ values were those obtained from the optimization of the previous time step. For the first excited state, 6 angles were optimized, resulting in an average of 7.45 iterations per nuclear time step over a single trajectory (with a time step of 0.5 fs). In the case of the second excited state, 13 angles were utilized, but only 7 were optimized, which required an average of 23.45 iterations. For the third excited state, 19 angles were involved, with 6 optimized, averaging 26.25 iterations per nuclear time step.

The distribution of coordinates and velocities was also generated by a Wigner distribution and the corresponding spectrum is shown in Fig. 7. The starting conditions (geometries, velocities, and initial state) were selected according to the excitation window ranging from 10.5 eV to 11.0 eV, which corresponds to the excitation to the S₂ state (ππ* excitation). A total of 100 initial conditions were extracted from a nuclear ensemble of 5000 geometries. In total, 85 were used to analyze the dynamics; 15 trajectories were discarded because the total energy was not conserved by the end of the dynamics. This energy inconsistency is primarily due to orbital rotation within the active space, resulting in discontinuities in the PESs and thus in violation in the total energy conservation. Capturing all different processes occurring in ethylene upon light irradiation, also including dissociation of hydrogen atoms, is a notorious-complicated problem that requires a large active space.⁸⁷ For the sake of practicality, we keep our small active space and removed the problematic trajectories. In the SH simulations, we used the local diabatization algorithm for the integration of the electronic TDSE, with an integration time step of 0.02 fs for the electronic degrees of freedom and 0.5 fs for the nuclear degrees of freedom. For some trajectories where the total energy was not conserved, the nuclear time step was reduced to 0.25 fs. The Granucci and Persico energy-based decoherence correction with the standard 0.1 a.u. parameter was applied during the SH dynamics. The rescaling of the nuclear velocities after a hop was performed in the direction of the nuclear momentum. The SH trajectories were propagated up to 100 fs; however, the trajectories that reached the ground state were stopped after running on the ground state for at least 10 fs.

Fig. 7 Simulated absorption spectrum of CH₂=CH₂ calculated using 5000 geometries. Contributions from each adiabatic state are indicated, with the dashed line representing the total spectrum.

3.2.2 Simulation of the excited states dynamics. The simulated absorption spectrum, shown in Fig. 7 exhibits one peak with two contributions. They come from the S₁ and S₂ states, even if they both correspond to the ππ* transition, which for different initial geometries is located at one of the S₁ or S₂ states. Most of the time, for the geometries picked by the Wigner distribution, the S₁ state corresponds to the σπ* excitation, which does not contribute to the spectrum. We note that the S₁ and S₂ states are very close in energy.

CH₂=CH₂ contains six atoms and therefore has twelve internal degrees of freedom that are all free to relax during the excited state dynamics simulations. Its compact structure, extremely fast dynamics (under 100 fs) through a conical intersection, and significant conformational flexibility have made this molecule a frequently used system for benchmarks and testing method developments, including mixed-quantum classical dynamics methods.^62–75,88

Previous SH-based studies of CH₂=CH₂ have discussed its dynamics in detail.^{62,63,69,70,73} In essence, after the ππ* state is excited, it decays very rapidly to the ground state. In doing so, CH₂=CH₂ transfers strong intramolecular vibrational redistribution of energy takes place, and besides photoisomerizing, it can show photodissociation. For example, it can release one or two H atoms and ethylidene isomer (CH₃–CH) is also a common photoproduct. Yet, CH₂=CH₂ is still the most common structure.⁶³ Several static studies and dynamics reveal that the most common conical intersection structure presents a twisted-pyramidalized geometry.^89–91

In our analysis, we do not focus on identifying and quantifying the photoproducts, as we terminate the dynamics after 10 fs in the electronic ground state.

In Fig. 8, we present two example trajectories. In the first trajectory (top), the S₁ state is populated just before 10 fs, and by 50 fs, the system undergoes a hop to the ground state. The hopping structure (geometry b1) exhibits a twisted-pyramidalized geometry, closely resembling the conical intersection geometry reported in previous static and dynamic studies.^{62–75,88–91} Consequently, the final geometry at the end of the dynamics (geometry c1) suggests that the final product is CH₂=CH₂.

Fig. 8 Potential energy surfaces calculated with VQE/VQD for two trajectories of CH₂=CH₂, illustrating (a1) and (a2) the initial configurations, (b1) and (b2) the hopping structure and (c1) and (c2) the final geometry. The thin black line below 2 eV represents the total energy (kinetic + potential) over time and the circles represent the active adiabatic state.

In the second trajectory (bottom of Fig. 8), we observe the decay of the S₂ state to the S₁ state shortly after 10 fs, followed by a rapid decay to the ground state just after 20 fs. In this case, the hopping structure (geometry b2) also displays a twisted-pyramidalized geometry. However, the final geometry appears to be in the process of forming the ethylidene isomer, as indicated by the apparent migration of a hydrogen atom to the other carbon.

Fig. 9 displays the population dynamics of the adiabatic states averaged over all trajectories following ππ* excitation. Initially, all populations reside in the S₂ state, characterized by ππ* character. The S₂ state undergoes a very fast exponential decay to the S₁ state which decays rapidly to the ground state. It is crucial to emphasize that in our dynamics, the σπ* state is never populated. Instead, the ππ* state, which initially corresponds to the S₂ and later becomes the S₁ state, decays to the ground state.

Fig. 9 Adiabatic state populations as functions of time, averaged over 85 trajectories for CH₂=CH₂. The gray line represents the fit of the S₁ and S₂ populations according to the kinetic model described by eqn (22).

Previous SH dynamics studies also reported that the initially populated ππ* state is responsible for decaying to the ground state.^62,63,69,70

To quantify the population decay, we fit a two-step irreversible kinetics model, defining the lifetimes τ₁ and τ₂ of the respective excited state manifolds as

where τ₂ is the lifetime of the S₂ population (P_S2), yielding τ₂ = 14 fs, and τ₁ is the lifetime of the S₁ population (P_S1), yielding τ₁ = 35 fs. This ultrafast decay of ∼50 fs (τ₁ + τ₂) aligns well with experimental results, where time-resolved spectroscopy has shown that ethylene's internal conversion occurs in only 10–30 fs.^92–96 We note that our obtained lifetime is also in agreement with other nonadiabatic dynamics studies performed for CH₂=CH₂ on classical computers.^62–75,88

4 Conclusions and future perspectives

We have reported a theoretical and computational framework to perform nonadiabatic dynamics simulations using a hybrid quantum-classical algorithm. By harnessing the variational quantum eigensolver and variational quantum deflation algorithms, we are able to accurately calculate ground and excited state energies, gradients, as well as nonadiabatic coupling vectors and transition dipole moments, which can all be used to perform nonadiabatic molecular dynamics simulations.

The trajectories were simulated within the framework of the Tully's fewest switches approach, using a local diabatization scheme for the integration of the electronic coefficients, albeit our method is versatile enough to be easily integrated into other mixed quantum-classical dynamics approaches. The presented methodology is integrated into our SHARC program package by making use of the TEQUILA quantum computing framework to compute the electronic properties.

Our approach was tested on two polyatomic molecular systems: the cis–trans photoisomerization of methanimine and the electronic relaxation of ethylene. The results demonstrate qualitatively accurate molecular dynamics for both cases, successfully reproducing findings from other computational studies and experimental data.

In this study, we employed the k-UpCCGSD ansatz for computing the electronic properties, but other ansätze can be seamlessly integrated. Furthermore, we anticipate extending this approach to other excited-state properties, such as spin–orbit couplings, by exploiting the flexibility of the VQD algorithm, which can include extra penalty terms to track spin multiplicity.

We believe this study represents a significant step towards realizing the “quantum advantage” for practical applications of quantum computers for chemical simulations. A thorough analysis of the robustness of the proposed algorithms against circuit noise present in real quantum computers is crucial before proceeding to run the calculations on NISQ quantum devices.

Data availability

All simulation results are included in the article and/or ESI.† The integration of the TEQUILA quantum computing framework into the SHARC program package will be available in an upcoming release of SHARC, which is published under the GPL license at https://github.com/sharc-md/.

Author contributions

E. S. G. conceptualized the study, implemented the methodology, performed dynamic simulations, curated the data, and wrote the original draft. M. O. supervised the project, assisted in the conceptualization, and contributed to the review and editing of the manuscript. J. S. K. assisted in the implementation and conceptualization, and contributed to the review and editing of the manuscript. L. G. supervised the project, acquired funding, and contributed to the review and editing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank the University of Vienna for generous continuous support. The computational results have been achieved in part using the Vienna Scientific Cluster (VSC).

Notes and references

1. R. P. Feynman, Int. J. Theor. Phys., 1980, 21, 467–488.
2. Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis and A. Aspuru-Guzik, Chem. Rev., 2019, 119, 10856–10915.
3. N. S. Blunt, J. Camps, O. Crawford, R. Izsák, S. Leontica, A. Mirani, A. E. Moylett, S. A. Scivier, C. Sünderhauf, P. Schopf, J. M. Taylor and N. Holzmann, J. Chem. Theory Comput., 2022, 18, 7001–7023.
4. R. Santagati, A. Aspuru-Guzik, R. Babbush, M. Degroote, L. González, E. Kyoseva, N. Moll, M. Oppel, R. M. Parrish, N. C. Rubin, M. Streif, C. S. Tautermann, H. Weiss, N. Wiebe and C. Utschig-Utschig, Nat. Phys., 2024, 20, 549–557.
5. P. Schleich, L. B. Kristensen, J. A. Campos-Gonzalez-Angulo, D. Avagliano, M. Bagherimehrab, A. Aldossary, C. Gorgulla, J. Fitzsimons and A. Aspuru-Guzik, arXiv, 2024, preprint, arXiv:2401.09268,  DOI:10.48550/arXiv.2401.09268.
6. J. Preskill, Quantum, 2018, 2, 79.
7. M. Brooks, Nature, 2019, 574, 19–21.
8. P. Adebayo, F. Basaky and E. Osaghae, UMYU Scientifica, 2023, 2, 65–75.
9. A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik and J. L. O'Brien, Nat. Commun., 2014, 5, 4213.
10. J. R. McClean, J. Romero, R. Babbush and A. Aspuru-Guzik, New J. Phys., 2016, 18, 023023.
11. J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth and J. Tennyson, Phys. Rep., 2022, 986, 1–128.
12. O. Higgott, D. Wang and S. Brierley, Quantum, 2019, 3, 156.
13. R. Carobene, S. Barison and A. Giachero, Quantum Sci. Technol., 2023, 8, 035014.
14. Q. X. Xie, S. Liu and Y. Zhao, J. Chem. Theory Comput., 2022, 18, 3737–3746.
15. J. Wen, D. Lv, M.-H. Yung and G.-L. Long, Quantum Engineering, 2021, 3, e80.
16. P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta and I. Tavernelli, Phys. Rev. Res., 2020, 2, 043140.
17. S. Zhao, D. Tang, X. Xiao, R. Wang, Q. Sun, Z. Chen, X. Cai, Z. Li, H. Yu and W.-H. Fang, J. Phys. Chem. Lett., 2024, 15, 7244–7253.
18. Y. Wang and D. A. Mazziotti, Phys. Chem. Chem. Phys., 2024, 26, 11491–11497.
19. I. O. Sokolov, P. K. Barkoutsos, L. Moeller, P. Suchsland, G. Mazzola and I. Tavernelli, Phys. Rev. Res., 2021, 3, 013125.
20. Q. Gong, Q. Man, J. Zhao, Y. Li, M. Dou, Q. Wang, Y.-C. Wu and G.-P. Guo, J. Chem. Phys., 2024, 160, 124103.
21. D. Hirai, W. Mizukami and Y. O. Nakagawa, Chem. Phys. Lett., 2023, 816, 140404.
22. A. Gandon, A. Baiardi, P. Ollitrault and I. Tavernelli, J. Chem. Theory Comput., 2024, 20, 5951–5963.
23. D. Bultrini and O. Vendrell, Commun. Phys., 2023, 6, 328.
24. S. Tamiya, S. Koh and Y. O. Nakagawa, Phys. Rev. Res., 2021, 3, 023244.
25. J. C. Tully and R. K. Preston, J. Chem. Phys., 1971, 55, 562–572.
26. S. Mai, D. Avagliano, M. Heindl, P. Marquetand, M. F. S. J. Menger, M. Oppel, F. Plasser, S. Polonius, M. Ruckenbauer, Y. Shu, D. G. Truhlar, L. Zhang, P. Zobel and L. González, SHARC3.0: Surface Hopping Including Arbitrary Couplings – Program Package for Non-Adiabatic Dynamics, 2023, https://sharc-md.org/.
27. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
28. J. S. Kottmann, S. Alperin-Lea, T. Tamayo-Mendoza, A. Cervera-Lierta, C. Lavigne, T.-C. Yen, V. Verteletskyi, P. Schleich, A. Anand, M. Degroote, S. Chaney, M. Kesibi, N. G. Curnow, B. Solo, G. Tsilimigkounakis, C. Zendejas-Morales, A. F. Izmaylov and A. Aspuru-Guzik, Quantum Sci. Technol., 2021, 6, 024009.
29. P. Jordan and E. Wigner, Z. Phys., 1928, 47, 631–651.
30. A. F. Izmaylov, T.-C. Yen, R. A. Lang and V. Verteletskyi, J. Chem. Theory Comput., 2019, 16, 190–195.
31. Z. P. Bansingh, T.-C. Yen, P. D. Johnson and A. F. Izmaylov, J. Phys. Chem. A, 2022, 126, 7007–7012.
32. T.-C. Yen and A. F. Izmaylov, PRX Quantum, 2021, 2, 040320.
33. A. Anand, P. Schleich, S. Alperin-Lea, P. W. K. Jensen, S. Sim, M. Díaz-Tinoco, J. S. Kottmann, M. Degroote, A. F. Izmaylov and A. Aspuru-Guzik, Chem. Soc. Rev., 2022, 51, 1659–1684.
34. F. A. Evangelista, J. Chem. Phys., 2011, 134, 224102.
35. I. G. Ryabinkin, T.-C. Yen, S. N. Genin and A. F. Izmaylov, J. Chem. Theory Comput., 2018, 14, 6317–6326.
36. J. Chen, H.-P. Cheng and J. K. Freericks, J. Chem. Theory Comput., 2022, 18, 2193–2198.
37. A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow and J. M. Gambetta, Nature, 2017, 549, 242–246.
38. D. A. Mazziotti, Phys. Rev. A: At., Mol., Opt. Phys., 2007, 75, 022505.
39. S. E. Smart and D. A. Mazziotti, Phys. Rev. Lett., 2021, 126, 070504.
40. S. Warren, Y. Wang, C. L. Benavides-Riveros and D. A. Mazziotti, Phys. Rev. Lett., 2024, 133, 080202.
41. J. Lee, W. J. Huggins, M. Head-Gordon and K. B. Whaley, J. Chem. Theory Comput., 2019, 15, 311–324.
42. J. J. Foley IV and D. A. Mazziotti, Mol. Phys., 2010, 108, 2543–2550.
43. Á. Ruiz-Serrano, N. D. M. Hine and C.-K. Skylaris, J. Chem. Phys., 2012, 136, 234101.
44. T. U. Helgaker and J. Almlöf, Int. J. Quantum Chem., 1984, 26, 275–291.
45. J. S. Kottmann, A. Anand and A. Aspuru-Guzik, Chem. Sci., 2021, 12, 3497–3508.
46. S. Gocho, H. Nakamura, S. Kanno, Q. Gao, T. Kobayashi, T. Inagaki and M. Hatanaka, npj Comput. Mater., 2023, 9, 13.
47. M. Barbatti, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 620–633.
48. R. Crespo-Otero and M. Barbatti, Chem. Rev., 2018, 118, 7026–7068.
49. S. Mai, P. Marquetand and L. González, in Surface Hopping Molecular Dynamics, John Wiley & Sons, Ltd, 2020, ch. 16, pp. 499–530.
50. L. Wang, J. Qiu, X. Bai and J. Xu, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2020, 10, e1435.
51. L. Verlet, Phys. Rev., 1967, 159, 98–103.
52. S. Fernandez-Alberti, A. E. Roitberg, T. Nelson and S. Tretiak, J. Chem. Phys., 2012, 137, 014512.
53. G. Granucci, M. Persico and A. Toniolo, J. Chem. Phys., 2001, 114, 10608–10615.
54. F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico and H. Lischka, J. Chem. Phys., 2012, 137, 22A514.
55. M. Barbatti and K. Sen, Int. J. Quantum Chem., 2016, 116, 762–771.
56. C. Pieroni, F. Becuzzi, L. Creatini, G. Granucci and M. Persico, J. Chem. Theory Comput., 2023, 19, 2430–2445.
57. Y. Yao, W. L. Hase, G. Granucci and M. Persico, J. Chem. Phys., 2021, 154, 074115.
58. L. Sun and W. L. Hase, J. Chem. Phys., 2010, 133, 044313.
59. D. Avagliano, E. Lorini and L. González, Philos. Trans. R. Soc., A, 2022, 380, 20200381.
60. Y. Suzuki, Y. Kawase, Y. Masumura, Y. Hiraga, M. Nakadai, J. Chen, K. M. Nakanishi, K. Mitarai, R. Imai, S. Tamiya, T. Yamamoto, T. Yan, T. Kawakubo, Y. O. Nakagawa, Y. Ibe, Y. Zhang, H. Yamashita, H. Yoshimura, A. Hayashi and K. Fujii, Quantum, 2021, 5, 559.
61. F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand and L. González, J. Chem. Theory Comput., 2016, 12, 1207–1219.
62. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys., 2020, 22, 15183–15196.
63. M. Barbatti, J. Chem. Theory Comput., 2021, 17, 3010–3018.
64. M. Ben-Nun, J. Quenneville and T. J. Martínez, J. Phys. Chem. A, 2000, 104, 5161–5175.
65. T. Mori, W. J. Glover, M. S. Schuurman and T. J. Martinez, J. Phys. Chem. A, 2012, 116, 2808–2818.
66. H. Tao, T. K. Allison, T. W. Wright, A. M. Stooke, C. Khurmi, J. van Tilborg, Y. T. Liu, R. W. Falcone, A. Belkacem and T. J. Martinez, J. Chem. Phys., 2011, 134, 244306.
67. J. K. Yu, C. Bannwarth, E. G. Hohenstein and T. J. Martínez, J. Chem. Theory Comput., 2020, 16, 5499–5511.
68. K. Saita and D. V. Shalashilin, J. Chem. Phys., 2012, 137, 22A506.
69. M. Barbatti, M. Ruckenbauer and H. Lischka, J. Chem. Phys., 2005, 122, 174307.
70. M. Barbatti, G. Granucci, M. Persico and H. Lischka, Chem. Phys. Lett., 2005, 401, 276–281.
71. M. Mališ and S. Luber, J. Chem. Theory Comput., 2020, 16, 4071–4086.
72. E. Pradhan, K. Sato and A. V. Akimov, J. Phys.: Condens. Matter, 2018, 30, 484002.
73. D. Hollas, L. Šištík, E. G. Hohenstein, T. J. Martínez and P. Slavíček, J. Chem. Theory Comput., 2018, 14, 339–350.
74. B. Sellner, M. Barbatti, T. Müller, W. Domcke and H. Lischka, Mol. Phys., 2013, 111, 2439–2450.
75. M. Fischer, J. Handt and R. Schmidt, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 90, 012527.
76. I. Tavernelli, U. F. Röhrig and U. Rothlisberger, Mol. Phys., 2005, 103, 963–981.
77. N. Minezawa and T. Nakajima, J. Chem. Phys., 2019, 150, 204120.
78. E. Fabiano, T. Keal and W. Thiel, Chem. Phys., 2008, 349, 334–347.
79. H. Hirai and O. Sugino, Phys. Chem. Chem. Phys., 2009, 11, 4570–4578.
80. N. L. Doltsinis and D. Marx, Phys. Rev. Lett., 2002, 88, 166402.
81. Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov and G. K.-L. Chan, J. Chem. Phys., 2020, 153, 024109.
82. G. Granucci and M. Persico, J. Chem. Phys., 2007, 126, 134114.
83. V. Bonačić-Koutecký and M. Persico, J. Am. Chem. Soc., 1983, 105, 3388–3395.
84. F. Schautz, F. Buda and C. Filippi, J. Chem. Phys., 2004, 121, 5836–5844.
85. A. Teslja, B. Nizamov and P. J. Dagdigian, J. Phys. Chem. A, 2004, 108, 4433–4439.
86. V. Bonačić-Koutecký and J. Michl, Theor. Chim. Acta, 1985, 68, 45–55.
87. L. Fransén, T. Tran, S. Nandi and M. Vacher, J. Phys. Chem. A, 2024, 128, 1457–1465.
88. W. Zhou, A. Mandal and P. Huo, J. Phys. Chem. Lett., 2019, 10, 7062–7070.
89. N. Minezawa and M. S. Gordon, J. Phys. Chem. A, 2009, 113, 12749–12753.
90. M. Ben-Nun and T. J. Martínez, Chem. Phys., 2000, 259, 237–248.
91. T. K. Allison, H. Tao, W. J. Glover, T. W. Wright, A. M. Stooke, C. Khurmi, J. van Tilborg, Y. Liu, R. W. Falcone, T. J. Martínez and A. Belkacem, J. Chem. Phys., 2012, 136, 124317.
92. P. Farmanara, O. Steinkellner, M. T. Wick, M. Wittmann, G. Korn, V. Stert and W. Radloff, J. Chem. Phys., 1999, 111, 6264–6270.
93. P. Farmanara, V. Stert and W. Radloff, Chem. Phys. Lett., 1998, 288, 518–522.
94. J. M. Mestdagh, J. P. Visticot, M. Elhanine and B. Soep, J. Chem. Phys., 2000, 113, 237–248.
95. V. Stert, H. Lippert, H. H. Ritze and W. Radloff, Chem. Phys. Lett., 2004, 388, 144–149.
96. K. Kosma, S. A. Trushin, W. Fuss and W. E. Schmid, J. Phys. Chem. A, 2008, 112, 7514–7529.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sc04987j

---