Implementation of quasiclassical mapping approaches for nonadiabatic molecular dynamics in the PySurf package

Abstract

Quasiclassical methods for nonadiabatic molecular dynamics, based on Mayer–Miller–Stock–Thoss mapping, are implemented in the open source computer package PySurf. This complements the implementation of surface hopping approaches performed in previous studies, and leads to a unified code that allows nonadiabatic dynamics simulations using various mapping approaches (Ehrenfest dynamics, the linearised semiclassical initial value representation, the Poisson-bracket mapping equation, the “unity” approach for the identity operator, the spin mapping, and the symmetrical quasiclassical windowing method) as well as different flavours of surface hopping (fewest-switches, Landau–Zener, and a mapping-inspired scheme). Furthermore, a plugin is developed to provide diabatic vibronic models as input in a sum-of-products form. This opens the way to the benchmark of different types of trajectory-based propagators on different models, against exact quantum dynamical simulations performed, e.g., by the multiconfigurational time-dependent Hartree method. Illustrative calculations, performed using the whole set of available propagators, are presented for different harmonic and anharmonic two-state models, exhibiting various degrees of correlation between vibrational modes.

---

1 Introduction

Nonadiabatic phenomena in photochemistry involve the dynamics of nuclear wave packets evolving in coupled potential energy surfaces (PES) associated to different electronic states. The computationally most efficient – and therefore most commonly used – approaches to simulate nonadiabatic molecular dynamics are methods where the nuclei follow independent, classical-like trajectories.^1–4

The main advantage of trajectory-based methods, in contrast to full quantum dynamical methods for nuclear motion,^5–7 is the fact that the electronic Hamiltonian does not need to be precalculated and fitted, but can be evaluated “on-the-fly” for varying nuclear geometries. A notable exception in quantum dynamics is the direct-dynamics variational multiconfigurational Gaussian (dd-vMCG) wave packet method,^8,9 that is based on a superposition of Gaussian wave packets that follow non-classical coupled trajectories. This approach allows on-the-fly nuclear quantum dynamics usually within a local harmonic approximation; however, due to strong nonlinearities in the equations of motion, the applications of vMCG have so far been limited to relatively small molecular systems.

On the other hand, the use of precalculated surfaces, such as linear or quadratic vibronic coupling models,^7,10,11 is not necessarily a prerogative of quantum dynamical methods. Analytical Hamiltonian models are widely used in trajectory-based simulations to extend the dynamical propagation to long time scales,^12,13 or to benchmark approximate methods against exact quantum mechanical results.^14–18

The methods explored in this work follow a general formalism in which only the electronic dynamics are formally treated with a rigorous quantum mechanical description, represented by a time-dependent electronic wavefunction

where the kets |n(Q)〉 indicate a finite basis of N electronic (adiabatic or diabatic) states, that in general depend on the nuclear coordinates Q. The complex coefficients C_n(t) are propagated via the time-dependent Schrödinger equation,where H^(QM)_nm(Q, P) are the elements of a quantum mechanical Hamiltonian matrix, that generally depend on nuclear coordinates Q = (Q₁, …, Q_F) and momenta P = (P₁, …, P_F).

The nuclear coordinates and momenta are initialised by sampling from an initial distribution and evolve according to independent Hamiltonian trajectories,

where M_κ is the mass associated to the κ-th degree of freedom and V_eff is an effective PES defining the forces responsible for the nuclear motion.

The use of an effective potential is an unavoidable consequence of the use of classical (i.e., independent) trajectories for the nuclear coordinates. The choice of V_eff is not unambiguous and has led to the formulation of a wide variety of methods based on independent trajectories.⁴ In general terms, the various approaches can be classified in two categories:

1. Mixed quantum-classical surface hopping methods, where V_eff(Q) coincides, at each time, with the PES of an active adiabatic state. Different approaches – such as the fewest switches,¹⁹ the Landau–Zener,²⁰ or the mapping-inspired^21,22 schemes – differ in the algorithm that determines the time step at which the active state should change and how the energy conservation should be imposed, i.e. when and how a trajectory should “hop” from one to another surface.

2. Quasiclassical mapping methods or “Ehrenfest-like approaches”, where the effective potential for the nuclear dynamics is obtained as a mean-field average of the electronic Hamiltonian. In contrast to surface hopping, the mapping methods can be rigorously derived as a classical limit of quantum mechanics,^23,24 where the real and imaginary part of the coefficients C_n(t) are proportional to “electronic” coordinates and momenta. Therefore, we denote these approaches as “quasiclassical” rather than “mixed quantum-classical”. The classical limit introduces an ambiguity in the sampling of the electronic phase space coordinates and in the evaluation of the electronic observables, leading to a variety of different mapping methods.^25–29

Notable recent developments on mapping methods include unified formulations for various approaches,^30,31 rigorous definitions of the zero-point energy parameter for the electronic variables,^31–33 as well as accurate windowing schemes to evaluate the electronic population.^34,35 Furthermore, ab initio implementations^36,37 and initial benchmark studies^15,17,38 have recently appeared.

Nevertheless, although the theoretical photochemistry community has gained significant expertise in surface hopping methods, mapping approaches have not been tested or applied as extensively. One reason is that using an average potential can seem nonphysical, particularly when nuclear trajectories move into regions where the potential energy surfaces of different states are well separated after passing through a near-degeneracy zone. Furthermore, these methods are known to break detailed balance in most cases,^39,40 a property necessary for internal consistency. The symmetric quasiclassical mapping approach by Cotton and Miller is partially an exception, since it has been shown to obey detailed balance in the limit of vanishing electron–nuclear coupling.⁴¹

Since mapping approaches constitute a proper classical limit of quantum mechanics, they are naturally more suitable for the development of non-classical corrections to impose microscopic reversibility,^39,40,42 including those based on interacting trajectories,^43,44 or to formulate quantum-classical analogues of stationary and time-resolved spectroscopic observables. In particular, trajectory-based approached for time-resolved nonlinear spectroscopy have been presented and reviewed by several authors.^45–48 Recent applications include simulations of the transient absorption and two-dimensional spectrum of pyrazine,^49,50 as well as pump–probe signals for azomethane⁵¹ and a dendrimer structure.⁵²

Another reason that limits the application of mapping approaches to nonadiabatic dynamics is the lack of a general-purpose code to systematically compare the various variants of different families of methods. Connected to this, a computational platform is missing to systematically benchmark surface hopping or quasiclassical mapping methods against numerically exact quantum dynamical results obtained, for example, by the Heidelberg MCTDH package^53,54 or Quantics.⁵⁵

Recently, the Pysurf package, developed by us, was presented.⁵⁶ This software is designed to facilitate prototyping and development tasks in general computational chemistry and, in particular, in the exploration of ground and excited state PES within nonadiabatic molecular dynamics. Due to its modular structure, it is an ideal platform to implement and test different propagation schemes for nonadiabatic dynamics.

The goal of this work is to extend the code to provide a computational platform for the systematic comparison between surface hopping and mapping methods, as well as to facilitate comparison with quantum dynamical calculations.

The implementation and testing of fewest-switches and Landau–Zener surface hopping schemes is described in a previous study.⁵⁷ Here we present the implementation of different propagation schemes based on the so called Mayer–Miller–Stock–Thoss mapping,^23,24 including a recently developed mapping-inspired approach to surface hopping.²² The code comes with a plugin to allow users to provide analytic PES in the sum of product form, using a format similar to that used for the MCTDH or Quantics codes. This enables the systematic and seamless benchmark of various mixed quantum-classical or quasiclassical methods against full quantum dynamical results.

The rest of the paper is organised as follows. Section 2 provides the essential derivation and an overview of the mapping methods for nonadiabatic dynamics implemented in this work, Section 3 describes the main details of the computational implementation, Section 4 illustrates prototypical applications where different methods are compared, and Section 5 summarises and concludes.

2 Mapping approaches for nonadiabatic molecular dynamics

Many variants of the nonadiabatic mapping approaches have been formulated over several decades, and derived using different formalisms.^{23–26,29,34,58,59} The purpose of this section is to provide a concise overview of the family of techniques implemented in the PySurf package, highlighting the fact they can all be connected to the classical limit obtained within the phase space formulation of quantum mechanics.

In the following we refer to nonadiabatic molecular dynamics that, in a full quantum mechanical setup, is described by the time-dependent Schrödinger equation

that accounts for an arbitrary number of electronic states (N) and nuclear degrees of freedom (F). In the diabatic representation, we take the molecular Hamiltonian in the formwhere |n〉 and |m〉 denote the discrete diabatic electronic levels, and the vectors Q = (Q₁,…,Q_F) and P = (P₁,…,P_F) collect the nuclear coordinates and momenta, respectively. V₀(Q) is an average diabatic PES and the matrix V(Q) = {V_nm(Q)} is defined to be real and traceless, i.e. .^28,30 For simplicity we define the initial state to be a nuclear wavefunction χ₀(Q) associated to one specific diabatic state, i.e.

|Ψ, t = 0〉 = |n〉χ₀(Q). (6)

The key step to derive quasiclassical mapping approaches for nonadiabatic molecular dynamics is the Schwinger mapping of an N-level system to the manifold of singly excited states of a N-dimensional harmonic oscillator,^24,60

where â_m and â^†_n are bosonic annihilation and creation operators and |0₁,…,1_n,…,0_N〉 is a state of a N-dimensional harmonic oscillator with a single excitation on the “electronic mode” n. Expressing the bosonic operators in terms of dimensionless coordinates and momenta, the Hamiltonian takes the so called Meyer–Miller–Stock–Thoss (MMST) mapping form,^23,24where q = (q₁,…,q_N) and p = (p₁,…,p_N) are the dimensionless coordinates and momenta for the mapping variables (in the Schrödinger representation p_n = −i∂/∂q_n). Combined with the time-dependent Schrödinger equationthe MMST Hamiltonian gives rise to exactly the same dynamics as eqn (4) and (6), provided that the initial state is defined consistently as a singly excited harmonic oscillator wavefunction

However, the quantum mechanical solution of eqn (9) with the initial condition given by eqn (10) is as difficult as solving the original problem. The potential advantage of the mapping approach lies in its straightforward formulation in the phase space formalism of quantum mechanics, upon which the classical limit can be defined rigorously. To this end, one introduces the Wigner function⁶¹

regarded as a quasi-probability distribution in the phase space, meaning that the expectation value of a generic observable Ω(Q, P, q, p) is evaluated asIn analogy with classical mechanics, the MMST Hamiltonian of eqn (8) is also treated as a function of the phase space variables, and not as an operator (therefore the caret is removed hereafter), and the time evolution of the Wigner function is given by the Wigner–Moyal equation,⁶²where the operator [capital Lambda, Greek, circumflex], consistently interposed between two phase space functions, is defined asand the arrow indicates the side on which the differentiation is done. Using the MMST Hamiltonian of eqn (8), the approximation e^{[capital Lambda, Greek, circumflex]} ≈ 1 + [capital Lambda, Greek, circumflex] gives the classical Liouville equation for the phase-space distribution,which can be solved by propagating classical (i.e., independent) trajectories according to Hamilton's equations:

2.1 Calculation of electronic observables

Eqn (16a)–(16d), whose propagation we implemented in PySurf, are common for most of the different mapping methods proposed over the last decades. The major differences between various approaches are the way the mapping variables q and p are sampled at the initial time and how the electronic populations are computed.

In our implementation the population in the state m, given that the dynamics is initiated in the state n, is evaluated using the so-called linearised semiclassical initial value representation (LSC-IVR),²⁵

where (q₀, p₀) and (q_t, p_t) are the mapping variables at the initial and final times and W_nuc(Q₀, P₀) is the initial distribution (Wigner function) for the nuclear coordinates. As usual, the integral of eqn (17) is evaluated by Monte Carlo sampling.

The functions F_nn and depend only on the mapping variables and define the phase space observables corresponding to the population operators |n〉〈n| and |m〉〈m|, respectively referred to the initial sampling and the evaluation of the electronic population at time t. As remarked by several researchers,^27,28 the use of the classical limit of eqn (15) makes the definition of such population functions ambiguous. Different choices for F_nn and , that would be formally equivalent in the full quantum treatment given by eqn (13), deliver different results when combined with classical trajectory dynamics. Indeed, one of the goals of the present implementation is to facilitate the comparison and the benchmark of different mapping approaches based on different definitions of F_nn and .

One possible choice is to map the electronic projector into a pure singly excited oscillator (SEO) state, and then evaluate the Wigner function of that state. This yields²⁸

whereNote that the quantity is conserved by the equations of motion (16c) and (16d).

Alternatively one can express the traceless part of the quantum mechanical projectors |n〉〈n| in terms of creation/annihilation operators, and then convert them into phase-space observables through a Wigner transformation,³⁰

In this case the sum over all the population functions is formally equivalent to the identity:^28,30 .

The functions F^SEO_nn(q, p) and F^Wig_nn(q, p) can be combined in all possible ways for the initial and final time in eqn (17), giving rise to different mapping methods, which have been proposed and derived independently in the literature. These approaches, implemented in PySurf, are listed below. For the following, it is useful to introduce the variables R_n² = q_n² + p_n², so that two different population functions can be expressed as

2.1.1 Linearised semiclassical initial value representation (LSC-IVR). In the original LSC-IVR approach²⁵ the Wigner transform of the singly excited harmonic oscillator wavefunction, given by eqn (20), is used both as initial distribution for the mapping variables and as a projector at time t. Formally, this implies that electronic populations are computed by replacing F_nn(q₀, p₀) = F^SEO_nn(q₀, p₀) and into eqn (17), where the normalisation factor (2π)^N results from the overlap between pure state Wigner functions.⁶² This giveswhere the last step follows from the fact the factor ϕ, being a function of R², is conserved during the dynamics (see eqn (19)). The last line of eqn (22) implies that the initial sampling over the mapping variables can be performed using the Gaussian distribution ϕ²(q(0), p(0)), and weighting the trajectories by the factor (2R_n²(0) − 1) in the evaluation of the observables. Although this is the most commonly used strategy, it has the disadvantage that the trajectories initially located near the circle R_n² = 1/2 get sampled frequently, but contribute to the ensemble average with a relatively low weight.

For this reason, in addition to the typical Gaussian sampling, we implemented the direct sampling from the absolute value Wigner distribution |(2R_n²(0) − 1)ϕ²(0)|. With this approach the trajectories are simply weighted by C × sign(2R_n²(0) − 1), where C is a normalisation constant. In the tests discussed in Section 4, this type of Wigner sampling leads to a faster convergence with respect to the number of trajectories.

2.1.2 Poisson bracket mapping equation (PBME). In the Poisson bracket mapping equation approach²⁶ the initial distribution for the mapping variables is given by F_nn(r₀, p₀) = F^SEO_nn(q₀,p₀), whereas the electronic populations are evaluated using the functions . With this choice, eqn (17) becomes

In this case the usual approach is to sample the initial trajectories from the Gaussian ϕ(q(0), p(0)) and evaluate the ensemble average by weighting the trajectories by the factor (2R_n²(0) − 1). To speed up the convergence, similarly to the LSC-IVR approach, we implemented the sampling from the absolute value distribution |(2R_n²(0) − 1)ϕ(0)|.

2.1.3 Improved population operator scheme. The schemes based on the so called “improved population operators”, formulated by Saller et al.,^28,30 are designed so that the total electronic population equals unity for each trajectory of the ensemble. One of such methods, denoted “single unity” in ref. 28 (or simply “unity” hereafter), is formally obtained by setting F_nn(q₀, p₀) = F^Wig_nn(q₀, p₀) and in eqn (17),where the last step follows from the conservation of ϕ. The equivalence between the expression of eqn (24) and that used by Saller et al.³⁰ is proved in the ESI.†

In our implementation, eqn (24) is evaluated by sampling from the Gaussian distribution ϕ(q(0), p(0)) and weighting the trajectories by the factor F^Wig_nn(0).

2.1.4 Spin mapping (SM). Another option to evaluate the population according to eqn (17) is to choose the F^Wig_nn functions for both the initial and final times. However, in the absence of the damping factor ϕ(q, p) the integral of eqn (17) diverges. This divergence can be eliminated by constraining the dynamics on the multi-dimensional spherical surface defined as , with a given fixed radius R̄ for all the trajectories. The justification for fixing the value of R² is the need to recover the Casimir invariant of the SU(N) group,⁵⁸ which describes quantum N-level systems.

To include the constraint, the mapping variables are sampled from the distribution F_nn(q, p) = AF^Wig_nn(q, p)δ(R² − R̄²), where the constant A = N!/(π^NR̄^2N−2) is fixed to normalise the distribution as . Evaluating the populations as

P^spin_m←n(t) = A〈F^Wig_mm(t)F^Wig_nn(0)δ(R²(0) − R̄²)W_nuc(0)〉 (25)

we can obtain the value of the sphere radius as ,⁵⁹ by simply requiring that at the initial time P_n←n(t) = 1. The derivation of the values of the constants A and R̄² is given in the ESI.†

The simulation of the electronic population dynamics using eqn (25) is equivalent to the spin mapping method introduced by Runeson and Richardson in the so called W-representation.^29,58 In their implementation the initial sampling is taken as uniform over the spherical surface. In PySurf the sampling is done directly from the absolute value distribution |F^Wig_nn(q(0), p(0))|δ(R²(0) − R̄²), and the trajectories are simply weighted by sign(F^Wig_nn(q(0), p(0))).

2.1.5 Ehrenfest. The standard Ehrenfest approach for nonadiabatic dynamics can be viewed as a special case of mapping approach, where the function F^Wig_nn is used for both the sampling at time t = 0 and the evaluation of the population at later times.

In a similar way to spin mapping, the integral of eqn (17) is made convergent by restricting the dynamics on a sphere. In this case the squared radius R̄² is taken to be equal to 2, by initialising the mapping variables taking R_n²(0) = 2 for the initial state n and R_m²(0) = 0 for the unpopulated states m. Formally, this corresponds to the expression

2.1.6 Symmetrical quasiclassical (SQC) mapping. The symmetrical quasiclassical windowing method developed by Cotton and Miller is also based on independent trajectories governed by the MMST Hamiltonian, but is formulated in terms of action-angle variables (see ref. 34 and 35 for details). In a nutshell, the approach is based on histogram-shaped population functions, that mimic the fact that the eigenvalues of the quantum mechanical observables R_n² are quantised and can take only odd integer values. To approximately recover this behaviour, Cotton and Miller proposed evaluating the populations using the expression

p^SQC_m←n(t) = A(t)〈F^SQC_mm(t)F^SQC_nn(0)W_nuc(0)〉, (27)

where histogram binning functions are used, i.e.and χ_A is the indicator function for the interval A. In our implementation, the parameter γ is set to , which is the recommended value for most applications.³⁵ In eqn (27) A(t) is a time-dependent normalisation constant, simply evaluated as³⁴In their study, Cotton and Miller proposed evaluating the electronic populations by using triangular binning functions, instead of histograms.³⁵ At present, these triangular functions have been used only in a few applications. Therefore, they will be implemented and tested in future work.

2.2 Adiabatic representation

Direct dynamics simulations, where the electronic Hamiltonian is obtained by quantum chemical calculations performed on-the-fly, are typically performed in the adiabatic representation, defined by the eigenstates of the electronic Hamiltonian.

Starting from the quantum mechanical MMST Hamiltonian of eqn (8), defined in the diabatic representation, the corresponding adiabatic Hamiltonian H̃_map is obtained by a unitary transformation

where is a unitary operator, appropriately chosen so that the electronic Hamiltonian matrix V(Q) is transformed to a diagonal form. To this end, we define the orthogonal matrix C(Q) of the coordinate-dependent eigenvectors of V(Q), and construct the adiabatic electronic Hamiltonian as

Ṽ(Q) = C^T(Q)V(Q)C(Q), (31)

where Ṽ(Q) is the diagonal matrix which contains the adiabatic PESs Ṽ_m(Q), which in the direct dynamics implementation are provided by electronic structure calculations. The adiabatic electronic states are obtained accordingly as linear combinations of diabatic states,

At this stage, we can define the diabatic-to-adiabatic unitary transformation in the quantum mechanical mapping formalism, as the transformation that changes the coordinates and momenta as

where the D^(κ)(Q) matrices are related to C(Q) byThese are anti-symmetric and contain the derivative couplings,

Note that the coordinate transformations of eqn (33a)–(33d) conserve the canonical commutators, which, according to Stone's theorem, guarantees the existence of the unitary operator which implements the transformation. With the above definition for , the adiabatic MMST Hamiltonian is obtained from eqn (8) and (30) as

In the classical limit, the nonadiabatic dynamics can be simulated by solving the related Hamilton's equations. To this end, it is convenient to propagate the kinematic momenta , instead of the canonical ones. With this choice the equations of motion take the formTherefore, dynamical simulations in the adiabatic representation do not require a diabatic-to-adiabatic transformation, but can instead be performed through direct computation of the adiabatic potentials Ṽ_n(Q) and the derivative couplings D^(κ)_nm(Q), i.e. the same quantities required for surface hopping simulations.

2.2.1 Mapping-inspired surface hopping. The mapping approaches described in the previous sections are mean-field approaches, where the nuclear dynamics is governed by an effective potential, obtained by averaging the diabatic or adiabatic PESs. This description has the well known disadvantage of becoming nonphysical when the nuclei move from a near-degeneracy zone towards regions where the potentials are largely separated in energy.

In contrast, surface hopping methods are based on propagating the nuclei on one “active” – typically adiabatic – PES. A possible strategy to define the active surface can be formulated from the adiabatic description given in Section 2.2, and leads to a so called mapping approach surface hopping (MASH). MASH is a non-stochastic technique based on jumps between active surfaces governed by the dynamics of the mapping variables q_n and p_n. The first formulation of this method was given by Mannouch and Richardson²¹ for two-state problems, and a multi-state “mapping-inspired” surface hopping (MISH) was later developed by Runeson and Manolopoulos for Hamiltonians with an arbitrary number of electronic states.²² Although a size-consistent multi-state generalisation of MASH has also been recently presented,⁶³ at present the implementation in PySurf follows the MISH formulation of ref. 22.

According to eqn (37b), the weight of each adiabatic potential to define the force on the nuclei is given by R_n²/2 = (q_n² + p_n²)/2. Starting from this observation, the key steps leading to MISH can be formulated as follows: (i) at each time during the dynamics, the active surface is defined as the one with the largest weight, so that the nuclear (kinematic) momenta evolve as²²

where(ii) the mapping variables are initialised in the sphere (so that the weights sum to one) and uniformly in the region where the prescribed initial state has the largest weight,with ;²² (iii) the populations are evaluated using a modified version of the F^Wig_mm function, where the traceless part is rescaled to obtain P_n←n(t = 0) = 1,To this end, α_N must be set equal to (N − 1)/(H_N − 1), where H_N is the N-th harmonic number.²² Note that with this definition it is guaranteed that the population functions sum to one.

As usual in surface hopping, whenever a change of active surface n→m occurs, the nuclear momenta need to be rescaled to ensure the conservation of energy. Following the derivation of Runeson and Manolopoulos,²² this is achieved by rescaling the momentum as

If, even with the rescaling, the conservation of energy cannot be ensured, the surface hop is aborted, and the momentum is reversed along the rescaling direction.

3 Implementation details

The architecture of the PySurf code, sketched in Fig. 1 was illustrated in detail in ref. 56. Briefly, the code is written in Python and consists of three parts. The core modules contain the Python classes that provide the basic functionalities for sampling, propagating classical trajectories, and getting PESs, forces and nonadiabatic couplings via the so called “surface-point-provider” (SPP).⁵⁶ The database modules include functionalities for storing and reading data (coordinates, energies, forces, etc.) using a specifically tailored engine based on the Network Common Data Form (NetCDF). Finally, the code comes with a Plugin engine that allows users to extend the code with specific sub-programs for sampling, propagating trajectories, interfacing electronic structure codes, interpolating precomputed data, or evaluating energies using predefined models.

Fig. 1 Architecture of the PySurf code: the core modules (blue) provide basic functionalities; the database modules (green) include classes to store and read coordinates, energies and dynamical data; the plugins (gray) allow adding new sub-programs with specialised tasks. The new plugins developed in this work are indicated in red. The version of PySurf used for the calculation of this work is freely available on GitHub.⁶⁴

In particular, the possibility of adding plugins is one of the key aspects that make the code modular and efficient.⁵⁶ The development of the present work involved the creation of three new sets of plugins, highlighted in red in Fig. 1: (i) samplers for the mapping variables, described in Section 2 for the various mapping methods; (ii) propagators to solve eqn (16a)–(16d), based on the integrator described in Section 3.1 below; (iii) a plugin to provide analytic potential energy surfaces in a general sum-of-product form, described in Section 3.2.

3.1 Integration of the quasiclassical mapping trajectories

Within the propagator plugins a new class implementing the coupled propagation of the nuclear and mapping variables was implemented. Given that the various approaches described in Section 2.1 rely on the same equations of motions, eqn (16a)–(16d), a unique integrator suffices for all of them.

The implemented integrator is inspired by the ones proposed in ref. 65. In the diabatic representation, for each time step, from t to t + δt, the following intermediate steps are performed:

1. The mapping variables are updated for a half-step (δt/2) by integrating eqn (16c) and (16d) while keeping the nuclear coordinates fixed. The integration is performed analytically as

2. The coordinates are integrated for the full step using the mapping variables at the half-step, i.e.,

3. The momenta are propagated for the full step using an average between the forces at the initial and final step, i.e.,

4. The mapping variables are integrated for the second half-step using the electronic Hamiltonian evaluated at time t + δt,

Similar equations are derived and implemented for the propagation in the adiabatic representation. Note that this scheme involves one evaluation of the electronic Hamiltonian and its gradients for each time step.

For computations based on surface hopping methods, the integration of the equations of motion is based on a velocity Verlet integrator, as described in previous studies.^57,66 In terms of computational cost, the propagation of a trajectory according to quasiclassical mapping methods is roughly as expensive as for surface hopping trajectories.

We noticed that for some mapping methods (in particular the PBME and the “unity” approaches) a relatively large number of trajectories are needed to converge the results. However, being classical, the trajectories can be propagated independently in parallel, which can be done in multi-node computing infrastructures, using the scripts provided with the code. On the other hand, running trajectories in parallel implies that the propagated observables (geometries, populations, etc.) are stored on disk for subsequent analysis. As described in ref. 56 this is done using the standard NetCDF3 database format, for which libraries are available in a number of programming languages. Unavoidably, a large number of trajectories and degrees of freedom correlates with a larger disk usage.

3.2 Input-given model electronic Hamiltonians

In order to facilitate the benchmark of the various mixed quantum-classical and quasiclassical methods against exact quantum mechanical results, the new class of SPP was implemented. This class provides the user with the possibility of defining analytical models, involving an arbitrary number of (diabatic) electronic states and potential energy surfaces and couplings, given as a sum-of-products (SOP) form.

The SOP form is crucial for the efficiency of full quantum dynamical methods based on tensor decomposition, such as those related to the multiconfigurational time-dependent Hartree (MCTDH) approach. With the new class the information about the electronic Hamiltonian governing the dynamics can be provided within a “potential energy surface file” that has a similar structure to the operator files used in the Heidelberg MCTDH package or in the Quantics code.^54,55 In this way the approximations derived from the classical treatment of the nuclear motion can be seamlessly verified, paving the way to systematic benchmarks for various classes of diabatic model, such as the linear or quadratic vibronic coupling models.^10,11

4 Example applications

This section illustrates the application of different trajectory-based approaches to various types of nonadiabatic processes, using predefined diabatic model Hamiltonians. In particular, simulations were performed using the six mapping approaches described in Section 2, as well as three different types of surface hopping methods, i.e., the Landau–Zener scheme (LZSH),²⁰ Tully's fewest-switches approach (FSSH),¹⁹ and the non-stochastic mapping-inspired surface hopping (MISH).²²

The goal of the present work is to illustrate the benefits of having various surface hopping and quasiclassical mapping methods implemented in the same code, which is instrumental to systematic benchmark studies. Such investigations should be preferably performed on larger sets of models for various classes of photochemical and photophysical processes,¹⁸ therefore they constitute the natural follow-up studies based on the present implementation.

Diabatic observables obtained by the trajectory-based calculations are compared with numerically exact full quantum results obtained by the MCTDH as implemented in the Heidelberg MCTDH package.⁵⁴ All simulations based on mapping approaches were performed in the diabatic representation. The surface hopping calculations were performed, as usual, in the adiabatic basis, and the computed observables were transformed to the diabatic basis afterwards, to appropriately compare them with the exact quantum wave packet results.¹⁴

The details of the computational setup of the calculations are provided in the ESI.† All input files and PySurf scripts used in this work are freely available in a GitHub repository.⁶⁴

4.1 B_2u → B_3u internal conversion in photoexcited pyrazine

The simulation of the internal conversion in pyrazine is one of the most common benchmark tests for validating new methods for nonadiabatic dynamics or new implementations. In this work we refer to the popular 4-modes linear vibronic coupling model that includes the excited diabatic states B_2u and B_3u. The quantum mechanical Hamiltonian is given aswith parameters taken from ref. 67. The trajectories are initialised in the diabatic state B_2u using the Wigner function , that corresponds to the ground vibrational state of the undisplaced harmonic oscillator.

The population of the B_3u state as a function of time, computed using different trajectory-based methods, is shown in Fig. 2. Panels (a) and (b) show the results obtained by various mapping approaches as well as the quantum mechanical reference (thick black line). In the numerically exact simulation, the B_3u states become populated up to ≈90% in 45 fs. At longer times the population exhibits recurrences where the population of the B_3u state reaches minima around 80 fs and 150 fs. These recurrence times are correctly predicted by all mapping approaches. However, all methods underestimate the extent of the short time (<50 fs) transfer. Nevertheless, the error is small (<10%) for the PBME, the SM and the SQC approaches, which anyhow recover the extent of population transfer quantitatively at longer times. The unity and Ehrenfest approaches are accurate at the recurrence times, but underestimate the maximum value of the B_3u population by more than 20%. The worst agreement with the reference result is obtained by the LSC-IVR that, although capturing the recurrence times, undervalues the extent of population transfer by 40% in the short time scale and by 20% at longer times.

Fig. 2 Time-dependent population of the B_3u state of pyrazine, calculated using exact quantum dynamics simulations (thick black line) and compared with (a) and (b) mapping approaches and (c) surface hopping approaches.

Fig. 2(c) compares the exact quantum mechanical results with the surface hopping predictions. All quantum-classical simulations predict the recurrence times correctly. In the first transfer event (<50 fs) all the surface hopping methods underestimate the B_3u population by ≈5% (for the FSSH and MISH methods) or by ≈10% (for the LZSH approach). At times longer than 70 fs the FSSH prediction agrees almost perfectly with the reference results, while both MISH and LZSH slightly underestimate the maximum population by 5–10%, with a MISH being decisively in better agreement with the reference. Errors are small at the recurrence times and larger when the B_3u population reaches its maximum.

A detailed investigation into the reasons why one method performs better than another requires extensive benchmarking across a large database of molecular models. Therefore, it is beyond the scope of the present work, whose main goal is to illustrate a novel implementation. Nevertheless, it is well known that mapping approaches – based on effective average potentials – often fail to describe wave packet branching, and this limitation is likely the source of the larger errors described above.

To illustrate this aspect, we plotted in Fig. 3 the vibrational distributions along the coupling mode Q₁ for the wave packet in the B_3u state, initially unpopulated. Panel (a) shows the numerically exact quantum mechanical result. The vibrational distribution is perfectly symmetric and exhibits a nodal line for Q₁ = 0, as a signature of the geometric phase effect. The density broadens and shrinks repeatedly in the range −5 < Q₁ < 5 with periods between 40 and 50 fs, branching in opposite directions with respect to the nodal line. This oscillatory dynamics is clearly discernible in the vibrational distributions computed using surface hopping methods [panels (h)–(j)], although the nodal line (a purely quantum feature) is not present and the neat oscillations slowly fade out at times >100 fs.

Fig. 3 Vibrational density distribution in the B_3u state as a function of the coupling mode Q₁ of pyrazine, computed (a) quantum mechanically, (b)–(g) using mapping approaches and (h)–(j) using surface hopping approaches. The atomic displacement vectors of the mode are sketched in the top right panel.

In contrast, the quasiclassical mapping methods [panels (b)–(g)] predict rather broad distributions with little structure. Clear oscillations are visible in the Ehrenfest results of panel (f). However, the amplitude of the oscillations is relatively small as compared to the exact quantum results. This is in line with the fact the Ehrenfest dynamics underestimates significantly the extent of B_2u → B_3u transfer. In comparison, a resemblance to oscillatory wavepacket branching is also noticeable in the SM and SQC results, which are the best performing mapping approaches.

Fig. 4 depicts the vibrational distribution along the tuning mode Q₂ – the most displaced mode in the B_3u surface – computed via quantum as well as trajectory-based methods. The quantum distribution, shown in panel (a), oscillates in the range −4 < Q₂ < 5 with a period of ≈60 fs, localising markedly at the turning points and exhibiting some wave packet branching. Also, for the dynamics of this mode, the best predictions are obtained by surface hopping methods [panels (h)–(j)]. An accurate description of the broadening and the localisation at the turning points is also obtained by the SQC and – to a lesser extent – the Ehrenfest methods [panels (g) and (f), respectively]. In contrast, the remaining mapping approaches significantly overestimate the amount of density around Q₂ ≈ 0. This qualitatively wrong result is in line with the worst performance in predicting the electronic population dynamics, especially as compared to the surface hopping methods.

Fig. 4 Vibrational density distribution in the B_3u state as a function of the tuning mode Q₂ of pyrazine, computed (a) quantum mechanically, (b)–(g) using mapping approaches and (h)–(j) using surface hopping approaches. The atomic displacement vectors of the mode are sketched in the top right panel.

Finally, it is worth noticing that for the LSC-IVR, unity and SM approaches, the vibrational distribution is locally negative (red areas in the plots). This is a well know problem of some mapping approaches, and it is due to the fact that the electronic mapping variables are sampled from a quantum mechanical (partially negative valued) Wigner function, which is then propagated using classical equations.

4.2 ²B₁ → ²B₂ internal conversion in the 2,6-bis(methylene) adamantyl cation

As an example of a high-dimensional system, we simulate the electronic relaxation from the ²B₁ to the ²B₂ states of the 2,6-bis(methylene) adamantyl (BMA) radical cation, modelled by the two-level full dimensional (78D) linear vibronic coupling model constructed in ref. 68. The diabatic model has the same form as that used for pyrazine, given in eqn (47), and was used in previous studies^17,38 to test various trajectory-based methods. The dynamics are initiated in the diabatic ²B₁ state using a Gaussian distribution for the vibrational coordinates, with the same form as that used in Section 4.1.

The population of the ²B₂ state as a function of time, computed using different surface hopping and mapping methods, is shown in Fig. 5. The reference quantum mechanical result was obtained by a multi-layer MCTDH calculation, whose details are given in the ESI.† In contrast to the pyrazine dynamics, the ²B₂ state is populated rather slowly (about 11% in 200 fs). Indeed, although the wave packet crosses a ²B₁/²B₂ intersection in 20–30 fs, it largely remains localised in a region where the energy gap between electronic states is relatively large (≈6400 cm⁻¹, see Fig. 10 of ref. 68 and the related discussion in the paper).

Fig. 5 Time-dependent population of the ²B₂ state of BMA, calculated using exact quantum dynamics simulations (thick black line) and compared with (a) and (b) mapping approaches and (c) surface hopping approaches.

As shown in panels (a) and (b), the population dynamics is correctly predicted by the PBME, unity and SM approaches until ≈150 fs. The LSC-IVR, unity and Ehrenfest methods are also relatively accurate, but overestimate the ²B₂ by ≈5% between 50 and 150 fs. At times >150 fs all mapping approaches deviate slightly from the exact quantum result. The present results are essentially identical to those of ref. 38, which were obtained using a different code; this reinforces the correctness of the present implementation.

The surface hopping results, shown in panel (c), compare differently with quantum dynamics, depending on the specific flavour of surface hopping. The FSSH and MISH approaches, that rely on the computation of nonadiabatic couplings, overestimate the population of the ²B₂ state by ≈5% throughout the entire simulation time, therefore they perform slightly worse than the best quasiclassical methods (PBME and unity). These results match nicely those of ref. 17 and 38.

The LZSH scheme has the worst performance. Landau–Zener trajectories, that are propagated only relying on adiabatic potential data, tend to hop completely to the lowest adiabatic surface at the first passage through the intersection (within 20–40 fs). Afterwards, no more transitions occur between adiabatic surfaces, and the diabatic populations start oscillating around the average value of ≈0.5 (not shown in the figure), as expected by the fact that the minima of the two diabatic surfaces differ only by 90 cm⁻¹.

Summarising, although more general conclusions require more extensive studies, the present results suggest that the quasiclassical mapping approaches might be more suitable to describe electronic relaxation dynamics in the case of weak diabatic coupling and rapidly separating surfaces.

4.3 B_3g→B_2g internal conversion in the butatriene cation

As a third example, we consider a two-state quadratic vibronic model for the coupled B_3g and B_2g state of the butratriene cation, that are accessible by photoionisation. In particular, we analyze the population dynamics resulting after initially populating the B_3g state with the same type of Gaussian Wigner function as for the case of pyrazine discussed in Section 4.1. For the simulation, the most general form of the diabatic quadratic vibronic coupling model is adopted,including all 18 modes of the molecule, with parameters taken from ref. 69.

The population dynamics predicted quantum mechanically and using various mapping and surface hopping approaches is depicted in Fig. 6. In the quantum calculation the B_2g state is populated up to 70% in the first 15 fs, and subsequently the population starts oscillating around the value of 0.6 with a period of ≈15 fs. The oscillations become more regular after 60 fs and their amplitude decays significantly after 120 fs. The long time population in the B_2g state is in the range 0.50–0.55.

Fig. 6 Time-dependent population of the B_2g state of the butatriene cation, calculated using exact quantum dynamics simulations (thick black line) and compared with (a) and (b) mapping approaches and (c) surface hopping approaches.

Among the mapping approaches, the PBME, the unity and the SM methods deliver very similar results. They all predict oscillations with the same period and amplitudes as for the quantum mechanical results, which persist also at a longer time scale. Although the population dynamics is reasonably accurate until 60 fs, these mapping approaches overestimate the B_2g population after 60 fs by ≈10%. Conversely, the LSC-IVR, Ehrenfest and SQC approaches underestimate the amplitude of the population oscillations, but predict the correct long time population. Interestingly, the LSC-IVR and Ehrenfest methods are the ones exhibiting the worst accuracy for the pyrazine model of Section 4.1.

In contrast to the mapping approaches, all surface hopping methods give similar results, as shown in Fig. 6(b). In particular, both FSSH and MISH are especially good in the first 60 fs. After that time, LZSH, FSSH and MISH deliver essentially identical predictions, slightly overestimating the regularity of the oscillations. Similar to the mapping approaches, LZSH and FSSH, they overestimate the long time B_2g population by a small extent (≈5%), whereas MISH is more accurate.

4.4 Photoisomerisation of a retinal chromophore model

The last example concerns the cis–trans photoisomerisation dynamics of the 2-cis-penta-2,4-dieniminium cation (cis-PSB3), a protonated Schiff base that can be regarded as a prototype model for the 11-cis retinal chromophore. PSB3 is described by the two-state three-coordinate anharmonic diabatic model constructed by Olivucci and coworkers.⁷⁰ The three coordinates are a bond length alternating vibration (r), a dihedral torsion (θ) that connects the cis and trans structures of PSB3, and a hydrogen out-of-plane wagging mode (ϕ).

Specifically, we use the diabatic model with the functional form and parameters given in ref. 16, neglecting the kinematic coupling between the various modes. This model is based on anharmonic functions that couple the three nuclear modes in a complex way, but nevertheless the total Hamiltonian has a SOP form that is suitable for MCTDH calculations. The plugin for arbitrary SOP diabatic potentials (see Section 3.2) allows running both quantum and quantum-classical simulations with essentially the same potential file. The initial condition for the nuclear coordinates is given by the same Gaussian distribution described in the Appendix of ref. 16, and the trajectories are initialised in the upper diabatic state (denoted state “2”) in a cis configuration.

In contrast to the previous examples, in this case the relevant quantity to be compared across different quantum-classical methods is the time-dependent quantum yield (TDQY) for the cis → trans isomerisation, computed as

where P⁽ⁿ⁾_cis/trans is the fraction of trajectories associated to a cis (cos θ > 0) or trans (cos θ < 0) configuration in the diabatic state n. Considering that, at time t = 0, P⁽²⁾_cis ≃ 1, the denominator of eqn (49) normalises over all possible photoproducts.¹⁶

Fig. 7 shows the TDQY computed quantum mechanically via MCTDH propagation, as well as using various quantum-classical approaches. The quantum dynamical simulation predicts a rise of the TDQY starting around 25 fs and peaking at ≈0.7 at 70 fs. This maximum is followed by a decay back to 0.4 at 120 fs. After this time, the TDQY oscillates back and stabilises around 0.5 with less pronounced oscillations. This behaviour differs from that of ref. 16, where a smaller primitive grid was used, but agrees nicely with the results of ref. 71.

Fig. 7 Time-dependent quantum yield for the formation of the trans-PSB3 isomer, calculated using exact quantum dynamics simulations (thick black line) and compared with (a) and (b) mapping approaches and (c) surface hopping approaches.

The quasiclassical mapping methods yield a different type of result. The LSC-IVR, PBME and unity approaches [panel (a)] fail at describing either the initial rise of the TDQY, and the subsequent decay and oscillatory stabilisation. Among them, the unity method is the only one able to predict the long time quantum yield quantitatively, whereas the LSC-IVR and PBME descriptions underestimate it by 5–10%. The SM, Ehrenfest and SQC methods [panel (b)] describe better the initial rise and decay of the TDQY, but fail in reproducing the oscillatory behaviour. This is likely due to the inadequacy of the mapping approaches in describing wave packet branching, as described for pyrazine in Section 4.1. Nevertheless, both the SM and SQC approaches predict the long time quantum yield with reasonable accuracy. Conversely, the Ehrenfest scheme overestimates the TDQY for times >50 fs.

Fig. 7(c) illustrates the predictions obtained by surface hopping. The three approaches, LZSH, FSSH and MISH, have a similar performance. FSSH replicates the TDQY curve obtained by quantum dynamics for the whole simulation time. The LZSH and MISH approaches respectively over- and under-estimate the extent of the initial rise in the quantum yield, but are remarkably accurate at times >100 fs, with MISH having a slightly better performance. The accuracy of the predictions for the diabatic population dynamics, shown in the ESI,† follows the same trends as for the TDQY.

It is worth noting that the error in the quantum yield for some approaches is relatively large compared to the ≈2% accuracy with which the quantum yield can nowadays be determined experimentally.⁷² This example shows that errors of ≈10% in the quantum yield, often attributed to an inadequate level of electronic structure theory, could also be due to a poor performance of certain trajectory-based methods for nonadiabatic dynamics. We hope that the present implementation will enable further benchmark studies to provide guidance on the choice of the “best” mixed quantum-classical or quasiclassical method to simulate molecular photoisomerisation dynamics.

5 Conclusions and perspectives

In this work various numerical methods based on the MMST mapping for nonadiabatic molecular dynamics are implemented. The implementation further includes the mapping-inspired surface hopping method by Runeson and Manolopoulos²² as well as a plugin to provide analytic diabatic potential energy surfaces in a sum-of-product form.

Building on previous implementations of the fewest-switches and Landau–Zener surface hopping schemes, the present development creates a comprehensive computational suite for quantum-classical nonadiabatic dynamics based on independent trajectories. All methods are seamlessly integrated into the Python package PySurf.⁶⁴ Its modular design enables the independent implementation of additional quantum-classical propagators, sampling algorithms, and evaluators for forces and nonadiabatic couplings, ensuring flexibility and ease of integration across different modules. In particular, users can easily implement other existing or novel quasiclassical approaches based on the MMST mapping. Following Fig. 1, this simply requires implementing a new sampling for the mapping variables, without the need for modifying either the propagator, the surface point provider or the interfaces with electronic structure codes.

In particular, the interface with different quantum chemistry packages, such as Q-Chem (used, e.g., in ref. 57), is automatically available in the new implementation of mapping propagators in the adiabatic representation (see Section 2.2). Such interface will be used in forthcoming investigations to run mapping and surface hopping calculations “on-the-fly”.

On the other hand, the plugin for analytic potential energy surfaces allows the use of essentially the same potential file to both PySurf and the most popular quantum dynamical codes, such as Heidelberg MCTDH or Quantics. This makes PySurf an excellent platform to perform systematic benchmarks of trajectory-based against full quantum dynamical results for large sets of Hamiltonian models and various classes of photochemical processes.¹⁸

The first tests for such benchmarks are presented here for prototypical photophysical and photochemical processes: the internal conversion in pyrazine, 2,6-bis(methylene) adamantyl and butatriene cation, as well as cis–trans isomerisation in a protonated Schiff base. Keeping in mind that the present simulations are far from constituting a comprehensive benchmark test, the current results suggest that the most commonly used FSSH method often provides the best agreement with quantum dynamics, together with the recently developed MISH scheme. Furthermore, in general, surface hopping schemes seem to require fewer trajectories for their convergence.

The perspective for future studies is to extend significantly the benchmark set, so to perform quantum and quantum-classical simulations systematically and on a high number of models. This will contribute to the creation of guidelines to help users to choose the best quasiclassical or quantum-classical approach to simulate specific photochemical phenomena.

Conflicts of interest

There are no conflicts to declare.

Data availability

All the computational results presented in this work can be reproduced using the data reported in the manuscript and the ESI.† The software used to perform the calculations, as well as the related input files, are freely available on GitHub: https://github.com/DavidPicconi/pysurf_dynamics/.

Acknowledgements

We would like to express our heartfelt gratitude to Prof. Christel M. Marian for her groundbreaking contributions to the advancement of theoretical and computational excited-state electronic structure methods. Her pioneering work, especially in the areas of multi-reference methods and spin-forbidden processes, has been instrumental in enhancing our understanding of light-responsive molecules. It is with great pleasure that we honour her remarkable achievements through the introduction of quasiclassical mapping approaches for nonadiabatic molecular dynamics in the PySurf package, featured in this special issue celebrating her 70th birthday. S. F. would like to express her gratitude to the Hans und Ria Messer Stiftung for their financial support of this work, under the project title ‘Hybride Quanten/Klassik Dynamiksimulationen photoaktiver Moleküle’ (project number Z_HMS20241).

References

1. M. Barbatti, Nonadiabatic dynamics with trajectory surface hopping method, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 620.
2. S. Mai, P. Marquetand and L. González, Nonadiabatic dynamics: the SHARC approach, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1370.
3. F. Agostini and B. F. E. Curchod, Different flavors of nonadiabatic molecular dynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1417.
4. S. Faraji, D. Picconi and E. Palacino-González, Advanced quantum and semiclassical methods for simulating photoinduced molecular dynamics and spectroscopy, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2024, 14, e1731.
5. Multidimensional Quantum Dynamics: MCTDH Theory and Applications, ed. H.-D. Meyer, F. Gatti and G. A. Worth, Wiley-VCH Verlag GmbH & Co. KGaA, 2009.
6. G. A. Worth, H.-D. Meyer, H. Koppel, L. S. Cederbaum and I. Burghardt, Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics, Int. Rev. Phys. Chem., 2008, 27, 569.
7. T. J. Penfold, E. Gindensperger, C. Daniel and C. M. Marian, Spin-Vibronic Mechanism for Intersystem Crossing, Chem. Rev., 2018, 118, 6975.
8. G. A. Worth, M. A. Robb and I. Burghardt, A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets, Faraday Discuss., 2004, 127, 307.
9. G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt and B. Lasorne, Quantum dynamics simulations using Gaussian wavepackets: the vMCG method, Int. Rev. Phys. Chem., 2015, 34, 269.
10. H. Koppel, W. Domcke and L. S. Cederbaum, Multimode Molecular Dynamics Beyond the Born-Oppenheimer Approximation, Adv. Chem. Phys., 1984, 59.
11. S. Faraji, S. Gomez-Carrasco and H. Koppel, Multistate Vibronic Dynamics and Multiple Conical Intersections, Conical Intersections: Theory, Computation and Experiment, World Scientific, 2011, p. 249.
12. F. Plasser, S. Gómez, M. F. S. J. Menger, S. Mai and L. González, Highly efficient surface hopping dynamics using a linear vibronic coupling model, Phys. Chem. Chem. Phys., 2019, 21, 57.
13. J. P. Zobel, M. Heindl, F. Plasser, S. Mai and L. González, Surface Hopping Dynamics on Vibronic Coupling Models, Acc. Chem. Res., 2021, 54, 3760.
14. W. Xie, M. Sapunar, N. Došlić, M. Sala and W. Domcke, Assessing the performance of trajectory surface hopping methods: ultrafast internal conversion in pyrazine, J. Chem. Phys., 2019, 150, 154119.
15. X. Gao, M. A. C. Saller, Y. Liu, A. Kelly, J. O. Richardson and E. Geva, Benchmarking Quasiclassical Mapping Hamiltonian Methods for Simulating Electronically Nonadiabatic Molecular Dynamics, J. Chem. Theory Comput., 2020, 16, 2883.
16. E. Marsili, M. Olivucci, D. Lauvergnat and F. Agostini, Quantum and Quantum-Classical Studies of the Photoisomerization of a Retinal Chromophore Model, J. Chem. Theory Comput., 2020, 16, 6032–6048.
17. Z. Liu, N. Lyu, Z. Hu, H. Zeng, V. S. Batista and X. Sun, Benchmarking various nonadiabatic semiclassical mapping dynamics methods with tensor-train thermo-field dynamics, J. Chem. Phys., 2024, 161, 024102.
18. L. E. Cigrang, B. F. E. Curchod, R. A. Ingle, A. Kelly, J. R. Mannouch, D. Accomasso, A. Alijah, M. Barbatti, W. Chebbi, N. Došlić, A. Freibert, L. González, G. Granucci, F. J. Hernández, J. Hernández-Rodríguez, A. Jain, J. Janoš, D. Lauvergnat, B. L. Dé, Y. Lee, S. K. Min, D. Picconi, E. S. Gil, M. Sapunar, P. Schürger, P. Vindel-Zandbergen, G. A. Worth, F. Agostini, S. Gómez, L. M. Ibele and A. Prlj, Roadmap for Molecular Benchmarks in Nonadiabatic Dynamics, J. Phys. Chem. A, 2025, 129(31), 7023–7050.
19. J. C. Tully, Molecular dynamics with electronic transitions, J. Chem. Phys., 1990, 93, 1061–1071.
20. A. K. Belyaev and O. V. Lebedev, Nonadiabatic nuclear dynamics of atomic collisions based on branching classical trajectories, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 014701.
21. J. R. Mannouch and J. O. Richardson, A mapping approach to surface hopping, J. Chem. Phys., 2023, 158, 104111.
22. J. E. Runeson and D. E. Manolopoulos, A multi-state mapping approach to surface hopping, J. Chem. Phys., 2023, 159, 094115.
23. H. Meyer and W. H. Miller, A classical analog for electronic degrees of freedom in nonadiabatic collision processes, J. Chem. Phys., 1979, 70, 3214–3223.
24. G. Stock and M. Thoss, Semiclassical Description of Nonadiabatic Quantum Dynamics, Phys. Rev. Lett., 1997, 78, 578–581.
25. X. Sun, H. Wang and W. H. Miller, Semiclassical theory of electronically nonadiabatic dynamics: results of a linearized approximation to the initial value representation, J. Chem. Phys., 1998, 109, 7064–7074.
26. H. Kim, A. Nassimi and R. Kapral, Quantum-classical Liouville dynamics in the mapping basis, J. Chem. Phys., 2008, 129, 084102.
27. A. Kelly, R. van Zon, J. Schofield and R. Kapral, Mapping quantum-classical Liouville equation: projectors and trajectories, J. Chem. Phys., 2012, 136, 084101.
28. M. A. C. Saller, A. Kelly and J. O. Richardson, On the identity of the identity operator in nonadiabatic linearized semiclassical dynamics, J. Chem. Phys., 2019, 150, 071101.
29. J. E. Runeson and J. O. Richardson, Spin-mapping approach for nonadiabatic molecular dynamics, J. Chem. Phys., 2019, 151, 044119.
30. M. A. C. Saller, A. Kelly and J. O. Richardson, Improved population operators for multistate nonadiabatic dynamics with the mixed quantum-classical mapping approach, Faraday Discuss., 2020, 221, 150–167.
31. J. Liu, X. He and B. Wu, Unified Formulation of Phase Space Mapping Approaches for Nonadiabatic Quantum Dynamics, Acc. Chem. Res., 2021, 54, 4215.
32. X. Gao and E. Geva, Improving the Accuracy of Quasiclassical Mapping Hamiltonian Methods by Treating the Window Function Width as an Adjustable Parameter, J. Phys. Chem. A, 2020, 124, 11006.
33. X. He, Z. Gong, B. Wu and J. Liu, Negative Zero-Point-Energy Parameter in the Meyer–Miller Mapping Model for Nonadiabatic Dynamics, J. Phys. Chem. Lett., 2021, 12, 2496.
34. S. J. Cotton and W. H. Miller, Symmetrical windowing for quantum states in quasiclassical trajectory simulations: application to electronically non-adiabatic processes, J. Chem. Phys., 2013, 139, 234112.
35. S. J. Cotton and W. H. Miller, A new symmetrical quasi-classical model for electronically non-adiabatic processes: application to the case of weak non-adiabatic coupling, J. Chem. Phys., 2016, 145, 144108.
36. D. Hu, Y. Xie, J. Peng and Z. Lan, On-the-Fly Symmetrical Quasi-Classical Dynamics with Meyer–Miller Mapping Hamiltonian for the Treatment of Nonadiabatic Dynamics at Conical Intersections, J. Chem. Theory Comput., 2021, 17, 3267.
37. B. M. Weight, A. Mandal, D. Hu and P. Huo, Ab initio spin-mapping non-adiabatic dynamics simulations of photochemistry, J. Chem. Phys., 2025, 162, 084105.
38. Y. Liu, X. Gao, Y. Lai, E. Mulvihill and E. Geva, Electronic Dynamics through Conical Intersections via Quasiclassical Mapping Hamiltonian Methods, J. Chem. Theory Comput., 2020, 16, 4479.
39. G. Amati, J. R. Mannouch and J. O. Richardson, Detailed balance in mixed quantum-classical mapping approaches, J. Chem. Phys., 2023, 159, 214114.
40. G. Amati, J. E. Runeson and J. O. Richardson, On detailed balance in nonadiabatic dynamics: from spin spheres to equilibrium ellipsoids, J. Chem. Phys., 2023, 158, 064113.
41. W. H. Miller and S. J. Cotton, Communication: note on detailed balance in symmetrical quasi-classical models for electronically non-adiabatic dynamics, J. Chem. Phys., 2015, 142, 131103.
42. H. W. Kim and Y. M. Rhee, Improving long time behavior of Poisson bracket mapping equation: a non-Hamiltonian approach, J. Chem. Phys., 2014, 140, 184106.
43. P. Huo and D. F. Coker, Communication: partial linearized density matrix dynamics for dissipative, non-adiabatic quantum evolution, J. Chem. Phys., 2011, 135, 201101.
44. S. K. Min, F. Agostini and E. K. U. Gross, Coupled-Trajectory Quantum-Classical Approach to Electronic Decoherence in Nonadiabatic Processes, Phys. Rev. Lett., 2015, 115, 073001.
45. I. Conti, G. Cerullo, A. Nenov and M. Garavelli, Ultrafast Spectroscopy of Photoactive Molecular Systems from First Principles: Where We Stand Today and Where We Are Going, J. Am. Chem. Soc., 2020, 142, 16117.
46. T. L. C. Jansen, Computational spectroscopy of complex systems, J. Chem. Phys., 2021, 155, 170901.
47. J. Krumland, M. Guerrini, A. De Sio, C. Lienau and C. Cocchi, Two-dimensional electronic spectroscopy from first principles, Appl. Phys. Rev., 2024, 11, 011305.
48. M. F. Gelin, Z. Lan, N. Došlić and W. Domcke, Computation of Time-Resolved Nonlinear Electronic Spectra From Classical Trajectories, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70012.
49. M. F. Gelin, X. Huang, W. Xie, L. Chen, N. Došlić and W. Domcke, Ab Initio Surface-Hopping Simulation of Femtosecond Transient-Absorption Pump–Probe Signals of Nonadiabatic Excited-State Dynamics Using the Doorway–Window Representation, J. Chem. Theory Comput., 2021, 17, 2394.
50. X. Huang, W. Xie, N. Došlić, M. F. Gelin and W. Domcke, Ab Initio Quasiclassical Simulation of Femtosecond Time-Resolved Two-Dimensional Electronic Spectra of Pyrazine, J. Phys. Chem. Lett., 2021, 12, 11736.
51. C. Xu, K. Lin, D. Hu, F. L. Gu, M. F. Gelin and Z. Lan, Ultrafast Internal Conversion Dynamics through the on-the-Fly Simulation of Transient Absorption Pump–Probe Spectra with Different Electronic Structure Methods, J. Phys. Chem. Lett., 2022, 13, 661.
52. D. Hu, J. Peng, L. Chen, M. F. Gelin and Z. Lan, Spectral Fingerprint of Excited-State Energy Transfer in Dendrimers through Polarization-Sensitive Transient-Absorption Pump–Probe Signals: On-the-Fly Nonadiabatic Dynamics Simulations, J. Phys. Chem. Lett., 2021, 12, 9710.
53. M. Beck, A. Jäckle, G. Worth and H.-D. Meyer, The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets, Phys. Rep., 2000, 324, 1.
54. G. A. Worth, M. H. Beck, A. Jäckle, H.-D. Meyer and O. Vendrell, MCTDH Package, Version 8.6, 2023, see https://mctdh.uni-hd.de.
55. G. Worth, Quantics: a general purpose package for Quantum molecular dynamics simulations, Comput. Phys. Commun., 2020, 248, 107040.
56. M. F. S. J. Menger, J. Ehrmaier and S. Faraji, PySurf: A Framework for Database Accelerated Direct Dynamics, J. Chem. Theory Comput., 2020, 16, 7681–7689.
57. E. X. Salazar, M. F. S. J. Menger and S. Faraji, Ultrafast Photoinduced Dynamics in 1,3-Cyclohexadiene: A Comparison of Trajectory Surface Hopping Schemes, J. Chem. Theory Comput., 2024, 20, 5796.
58. J. E. Runeson and J. O. Richardson, Generalized spin mapping for quantum-classical dynamics, J. Chem. Phys., 2020, 152, 084110.
59. D. Bossion, W. Ying, S. N. Chowdhury and P. Huo, Non-adiabatic mapping dynamics in the phase space of the SU(N) Lie group, J. Chem. Phys., 2022, 157, 084105.
60. J. Schwinger, in Quantum Theory of Angular Momentum, ed. L. C. Biedenharn and H. V. Dam, Academic, New York, 1965.
61. In eqn (11) and (14) it is implied that rectilinear coordinates are used for the nuclei, and the commutator between positions and momenta is [Q_κ, P_κ′] = iħδ_κκ′. In general, for rectilinear coordinates, the term ħ appearing in eqn (11) and (14) should be replaced by |[Q_κ, P_κ]|.
62. T. Curtright, D. Fairlie and C. Zachos, Features of time-independent Wigner functions, Phys. Rev. D: Part. Fields, 1998, 58, 025002.
63. J. E. Lawrence, J. R. Mannouch and J. O. Richardson, A size-consistent multi-state mapping approach to surface hopping, J. Chem. Phys., 2024, 160, 244112.
64. PySurf: A Framework for Database Accelerated Quantum Chemistry. See https://github.com/DavidPicconi/pysurf_dynamics.
65. J. O. Richardson, P. Meyer, M.-O. Pleinert and M. Thoss, An analysis of nonadiabatic ring-polymer molecular dynamics and its application to vibronic spectra, Chem. Phys., 2017, 482, 124.
66. E. Salazar Quezada, Mixed quantum/classical dynamics simulations of molecular excited processes, PhD thesis, University of Groningen, 2023.
67. G. A. Worth, H.-D. Meyer and L. S. Cederbaum, Relaxation of a system with a conical intersection coupled to a bath: a benchmark 24-dimensional wave packet study treating the environment explicitly, J. Chem. Phys., 1998, 109, 3518.
68. A. F. Izmaylov, D. Mendive–Tapia, M. J. Bearpark, M. A. Robb, J. C. Tully and M. J. Frisch, Nonequilibrium Fermi golden rule for electronic transitions through conical intersections, J. Chem. Phys., 2011, 135, 234106.
69. C. Cattarius, G. A. Worth, H.-D. Meyer and L. S. Cederbaum, All mode dynamics at the conical intersection of an octa-atomic molecule: multi-configuration time-dependent Hartree (MCTDH) investigation on the butatriene cation, J. Chem. Phys., 2001, 115, 2088.
70. E. Marsili, M. H. Farag, X. Yang, L. De Vico and M. Olivucci, Two-State, Three-Mode Parametrization of the Force Field of a Retinal Chromophore Model, J. Phys. Chem. A, 2019, 123, 1710–1719.
71. M. Aarabi, E. Marsili, M. Olivucci, D. Lauvergnat, F. Agostini, M. Garavelli and F. Santoro, Quantum Dynamics Predicts Coherent Oscillations in the Early Times of a Biological Photoisomerization, J. Phys. Chem. Lett., 2025, 16, 8486–8494.
72. M. Gueye, M. Paolino, E. Gindensperger, S. Haacke, M. Olivucci and J. Léonard, Vibrational coherence and quantum yield of retinal-chromophore-inspired molecular switches, Faraday Discuss., 2020, 221, 299.

---

Footnote

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

---