Direct wavepacket dynamics with spin–orbit coupling: simulation of thioformaldehyde

Abstract

The non-radiative photophysical phenomenon known as intersystem crossing is simulated, for the first time, in fully coupled direct nuclear wavepacket dynamics based on Gaussian basis functions (DD-vMCG method). As a case-study, the ultrafast photoinduced dynamics of thioformaldehyde is studied using this method. Working in the spin-diabatic basis, the simulations presented here incorporate the effect of spin–orbit coupling, thus allowing the transfer of population density between states of different spin multiplicities (in this case, singlets and triplets). Even though no significant intersystem crossing is expected in the first ps after exciting the first singlet excited state of thioformaldehyde, our calculations show that the subtle effects of non-negligible spin–orbit coupling elements are accounted for. Additionally, internal conversion back to the ground state is observed, mediated by non-adiabatic effects. While our reference simulations are based on MS-CASPT2 and MRCI calculations, we also investigate the sensitivity of the dynamics to different choices of underlying electronic structure methods. Finally, a comparison is made between results presented here, and those obtained from previously published surface-hopping simulations, commenting on possible methodological reasons for any discrepancies.

---

1 Introduction

In the fields of photophysics and photochemistry, the phenomenon of intersystem crossing (ISC) is well documented.^1–3 It involves the non-radiative transition of population between two electronic states of different spin multiplicities, making it analogous to internal conversion (IC) in which transitions between states of the same spin take place. Formally, a change in spin quantum number is forbidden by selection rules, but ISC was nevertheless identified to occur in certain molecules by Kasha, and others.^4,5 It is now known that the selection rule can be broken by an interaction between the angular momentum of electrons in orbitals, and their spin. This interaction is a relativistic effect referred to as spin–orbit coupling (SOC), and is of great academic and commercial interest. Specifically, the transitions from singlet to triplet states have been widely studied in the development of functional chromophores, photosensitisers, and optoelectronics.^6–8

While the strength of SOC is enhanced by heavier elements, the effect has also been documented in lighter small organic molecules, with a sulphur substituent for example.^9,10 One such molecule is thioformaldehyde (CH₂S, i.e. sulphur-substituted formaldehyde, hereafter referred to as TF), which has been used in several previous computational investigations to study the mechanism behind singlet–triplet interconversion.^11–17 Its small size makes it a perfect candidate for high-level computational studies, and hence, in this paper it will be used to show how SOC and its effect can be included in a fully quantum mechanical treatment of the nuclear dynamics of molecules.

In our simulations, we include four electronic states of TF; the ground state, the first singlet excited state (¹nπ*), and the first and second triplet states (³nπ* and ³ππ* respectively). The molecule will be excited from the ground state to the excited singlet state, after which the wavepacket created on that state is evolved in time using a quantum dynamical algorithm.

Quantum dynamics (QD) methods simulate molecules by propagating their wavefunction according to the time-dependent Schrödinger equation (TDSE), allowing for an accurate description of, for example, excited electronic state processes. By not relying on the Born–Oppenheimer approximation, the coupled motion of nuclear and electronic degrees of freedom (DOF) can be computed and non-adiabatic effects observed. Such simulations are the most complete way of describing excited state dynamics and hence are perfectly suited to study ISC. However, this accuracy comes at a great computational cost, scaling exponentially with system size. Naturally, this fact has led to the development of many approaches to solve the TDSE, some more approximate than others, and each with different characteristic features. In this work we describe one such method, whose origins can be traced back to the well-known and very powerful multi-configurational time-dependent Hartree (MCTDH) method.^18–21

In MCTDH, the wavefunction (and Hamiltonian) is represented by a number of nuclear basis functions that populate a grid. Such an implementation makes the method numerically efficient and able to converge on the exact solution of the TDSE. Its main drawback, common to all grid-based QD methods, is that a global description of the underlying potential energy surface (PES) is required, which can be prohibitively expensive. The procedure employed in this work does away with this requirement by calculating the potential only in those regions of configuration space that are actually visited by the moving basis functions. This is called direct (or on-the-fly) dynamics and is achieved, in this case, by expressing the nuclear basis set using Gaussian functions, hence its name; direct dynamics variational multi-configurational Gaussian (DD-vMCG).^22–25

The option of explicitly including the effects of SOC has only recently become available for DD-vMCG, implemented in the Quantics package.²⁶ The addition of this feature is an important step towards (1) making Quantics a more complete package for the study of quantum dynamical processes, and (2) allowing the DD-vMCG method to be compared to other QD methods with respect to describing ISC phenomena. To expand on this second point, we refer to a community-wide effort to produce a roadmap towards benchmarking non-adiabatic dynamics methods, which was recently published as a perspective article.²⁷ In it, four key photochemical/photophysical phenomena were identified to be of central interest, one of which being non-reactive radiationless relaxation (NRR). This includes ISC and thus illustrates the importance of being able to simulate it, in addition to its more commercial interest in the development of functional materials, as mentioned. IC is of course also a pathway to NRR, and it will therefore also be discussed in the context of TF. We note that IC is often ignored in studies on TF since the focus lies on ISC, but it is generally important to include all possible relaxation channels in the case of competing mechanisms.

2 Nuclear dynamics: the vMCG method

Following quantum theory, a complete description of molecular motion is given by the time-dependent Schrödinger equation (TDSE). Our purpose in this work is to describe spin–orbit coupling effects in molecular systems which implies the possibility of ISC, an inherently non-adiabatic phenomenon, and thus the need for a method which solves the TDSE while incorporating the coupling between nuclei and electrons. The variational multi-configurational Gaussian (vMCG) method meets those requirements by employing a coupled set of N time-dependent Gaussian basis functions, G, which describe a molecular wavepacket:where q are the nuclear coordinates, here taken to be the set of normal modes {q_k}, and N_s is the number of electronic states included. The multi-dimensional basis functions, referred to as Gaussian wavepackets (GWP), consist of a product of separable one-dimensional Gaussianswhich are each expressed in terms of a width, α, a linear parameter, ξ, dictating the position and momentum of the centre of the wavepacket, and a scalar parameter, η, containing the phase:

g^(k)_i(q_k,t) = exp(−α^(k)_iq_k² + ξ^(k)_i(t)q_k +η^(k)_i(t)). (3)

Notice that the widths are time-independent, implying the frozen Gaussian approximation. Additionally, by keeping the functions normalised, the scalar parameters no longer carry the phase information (which is instead given by the expansion coefficients of eqn (1)), and hence the time-dependence of the GWPs is carried entirely by the linear parameters.

The equations of motion (EoM) for the vMCG wavefunction ansatz are derived by applying the Dirac–Frenkel variational principle to both the coefficients and the Gaussian parameters, leading to coupled EoMs for each, expressed here in compact matrix notation:

iȦ = S⁻¹(H − iτ)A, (4)

and

i[capital Lambda, Greek, dot above] = C⁻¹Y. (5)

The coefficients collected in A are controlled by the overlap matrix, S_ij(t) = 〈G_i(t)|G_j(t)〉, differential overlap matrix, , and Hamiltonian matrix, H_ij = 〈G_i(t)|Ĥ|G_j(t)〉. The Hamiltonian operator, Ĥ, is expressed as a sum of the nuclear kinetic energy term, , and the electronic part, Ĥ^el, which defines the electronic states of interest and any interaction between them. The set of parameters for the Gaussian centres, Λ_i = {α_i,ξ_i,η_i}, are propagated according to two matrices, C and Y, which will not be discussed in detail. It will simply be noted that both contain the density matrix, , and the Gaussian overlap matrix. The former ensures the variational propagation of the Gaussian functions, while the latter is required because of the non-orthogonality of the basis set. The key point is that the GWP basis functions follow variationally coupled trajectories. All additional detail can be found in the original literature.^22,24

2.1 Interstate coupling

We now move to an overview of how the coupling between pairs of electronic states is included in vMCG. As a starting point, a distinction is made between same-spin coupling and coupling between states with different spin multiplicities. As such, the total electronic Hamiltonian, Ĥ^el, is split up in a non-relativistic (NR) part to account for the former, and a SOC contribution for the latter

Ĥ^el = Ĥ^NR + Ĥ^SOC. (6)

The second term is often ignored in simulations of small organic molecules where typically only singlet states are considered. In the adiabatic representation, the non-relativistic Hamiltonian matrix contains the potential energies of the electronic states (given by the eigenvalues, V^adia, of the spin pure electronic eigenstates, |ϕ〉) as on-diagonal elements, with coupling between states due to the derivative coupling terms, F, which are part of the kinetic energy operator (referred to as non-adiabatic coupling, NAC)²⁸

F_ss′ = 〈ϕ_s|∇ϕ_s′〉. (8)

Unfortunately the adiabatic representation cannot be used in vMCG simulations due to the singularities encountered in the derivative couplings at (near) degenerate points on the PESs (e.g. conical intersections).^29–31 To avoid these, the states require a unitary transformation to a diabatic basis, in which the couplings are smooth potential-like terms.²⁸ These states however, are now no longer eigenstates of the electronic Hamiltonian. The specific diabatisation scheme used here is discussed in the next section.

The main aim of this paper is to present, for the first time, the inclusion of the relativistic spin–orbit coupling Hamiltonian in DD-vMCG simulations. The SOC Hamiltonian matrix has the following elements, where the spin multiplicity of a state s is given by the subscript “m_s” = 2S + 1:

This additional coupling term can be explicitly added in the nuclear Schrödinger equation in two distinct ways, giving rise to spin diabatic and spin adiabatic representations.^32,33 The spin diabatic representation conserves the spin purity of electronic states and simply introduces the SOC as off-diagonal elements in the system Hamiltonian matrix. The Hamiltonian thus can be arranged with blocks of states of a single spin. For example if a system contains singlet (S) and triplet (T) states, then

with both singlet and triplet states rotated into a diabatic representation. In contrast, the spin adiabatic representation folds the SOC elements into a diagonalised Hamiltonian, giving rise to potentials for mixed spin states with coupling between them due only to the derivative couplings. In the following we employ the spin diabatic basis, which is both numerically more straightforward and chemically more intuitive. It is worth noting however that the use of one spin representation over another is not a straightforward choice, and much debate is still going on about their relative applicability in different regimes. We provide some general comments on this in the next section and in the SI, but unfortunately, direct comparison of spin diabatic versus spin adiabatic DD-vMCG calculations are impossible at the present time since the latter is not currently implemented.

2.2 On-the-fly propagation

Here we describe the direct dynamics (DD) extension to vMCG, used for all simulations presented in this work. This DD-vMCG method falls in the category of on-the-fly nuclear dynamics methods since the PES is gradually built up during the propagation of the wavepackets.²⁴ The need for pre-defined, fitted surfaces required for grid-based QD methods is thus avoided. The PES is locally expanded around the centre of each GWP, q₀, according to a second-order Taylor expansion known as the local harmonic approximation (LHA):

This expression allows for the analytical evaluation of the matrix elements of eqn (4) and (5) using information readily obtained from widely accessible quantum chemistry (QC) calculations.

Starting with this local description as a reference, new QC points will be calculated and added to a database as the GWPs move away from any points already in the database. A Shepard interpolation is then used to provide the LHA between the complete set of points in a database, essentially yielding “global” PESs that get more and more accurate as the number of records in the QC database increases.³⁴ For this reason it is practically advisable to perform multiple sequential DD-vMCG simulations for a given system, starting each time from the previously generated database. This ensures superior initial conditions for each additional run since the propagation is started on increasingly accurate surfaces.

Whenever a new point is required for the database, a quantum chemistry calculation is used to provide the adiabatic energies, V, of the excited states included in the simulation along with the nuclear gradients, V′, the non-adiabatic couplings between states of the same symmetry, F_ij(V_j − V_i), and the spin–orbit couplings between states of different spin, H^SOC. The Hessian matrices, V″, are explicitly calculated at the first point only and the curvature at other data points provided by Hessian updating.³⁵

The adiabatic data obtained is then transformed to a diabatic representation, used during the dynamics. The propagation diabatisation scheme³⁶ is employed in which the adiabatic-to-diabatic transformation matrix, D, is propagated along with the basis functions. An updated transformation matrix can be defined at any point q + Δq by the following relationship³⁷

∇D = −FD, (12)

using a previously obtained matrix D(q) and the derivative coupling matrix F_ss′. Such a scheme obviously requires an initial definition of D, for which we take the ground-state equilibrium geometry (i.e. the Franck–Condon, FC, point). This point serves as a reference (or gauge) where the adiabatic and diabatic representations are defined to be exactly equal. We note that when the number of electronic states is truncated this scheme does not produce an exact diabatisation (sometimes referred to as “quasi-diabatic”), but it does provide a very straightforward way to represent any number of adiabatic states in the required diabatic form. The generated PESs stored in the database can then be inspected to determine the quality of the diabatisation, and any necessary refinements can be made. In the present case, all electronic state surfaces for TF are successfully diabatised using the approached outlined above, without any adjustments.

The spin–orbit couplings are also rotated into the diabatic representation. For a full calculation the triplet states should be represented as a set of three sub-states with complex spin–orbit couplings. However it is found in practice that calculated spin–orbit coupling surfaces are often not smooth due to the random phase of the wavefunctions in quantum chemistry calculations³⁸ which makes them hard to use directly in simulations. In this work we therefore represent a triplet state in the commonly used form of a single averaged state, with real spin–orbit couplings given by the magnitude of the components:

While this simple effective SOC approach is known to have potential problems³³ and must be used with care, in particular if the SOCs show strong geometry dependence, or if more than one singlet state couples to the triplet state. We note that the nature of interstate couplings in TF, however, poses no fundamental problems for the chosen spin representation as these complications are not present. Furthermore, the validity of a spin representation may be sensitive to the chosen dynamics method, i.e., how the nuclei and possible state transitions are described,³⁹ and how these differences in representation transpire for DD-vMCG remains an open question. A more detailed discussion is provided in Section S5 of the SI.

3 Computational details

3.1 Electronic structure

One of the most crucial steps in any molecular dynamics study is the choice of electronic structure method with which the PESs will be calculated, and this is especially true for direct dynamics methods. In this case, we employ two closely related families of methods; multi-configurational self-consistent field (MCSCF) and multi-reference configuration interaction (MRCI). The former comprises the well-known (state-averaged) complete active space self consistent field (CASSCF) method^40,41 and its second-order perturbation theory extension, CASPT2.^42,43 A multi-state CASPT2 calculation serves as one of the “reference” calculations in this work and includes 4 electronic states (2 singlets, 2 triplets), calculated using 12 electrons in 10 active orbitals and employing the aug-cc-pVDZ basis set. A slightly modified active space, in which one occupied orbital was removed and one unoccupied orbital added, was employed in a MRCI(10,10)⁴⁴ calculation, using an ano-rcc-pVDZ basis set, and serves as a second reference (see the SI for the shape/descriptions of the orbitals).

In addition to these two references, various other CASSCF, CASPT2 and MRCI calculations were run to investigate how sensitive the results of our simulations are to the PESs calculated at different levels of theory. This sensitivity has been the topic of several recent works and is arguably one of the main bottlenecks in the field.^15,45–50 (as discussed in Section 2.3 of the benchmarking roadmap²⁷). Ultimately though, the effects of different electronic structure methods will be highly system specific – and to some extent, sensitive to the chosen dynamics method. A table with all electronic structure methods used in direct dynamics simulations is given in the SI, along with the predicted excitation energies for the first singlet (S₁), and two lowest triplet (T₁, T₂) states.

3.2 Direct dynamics

All QD calculations are carried out using the DD-vMCG method implemented in the Quantics package.^26,51 The coordinates used are the mass-frequency scaled ground-state normal modes of TF, calculated using an appropriate level of electronic structure theory, details of which are given in Table 1. A more elaborate comparison follows in Table S4 of the SI.

Table 1 Normal mode frequencies (in cm⁻¹) of TF calculated using both CCSD (Gaussian) and MRCI(10,10) (Molpro); utilised in reference calculations – MS-CASPT2(12,10) (Molcas) and MRCI(10,10) (Molpro), respectively

Label	CCSD	MRCI(10,10)	Description
q1	1028.09	969.52	CH2 oop bend
q2	1044.82	1002.10	CH2 rock
q3	1100.30	1021.26	C–S stretch
q4	1557.24	1484.21	CH2 scissor
q5	3134.95	3002.59	CH2 sym. stretch
q6	3222.79	3094.51	CH2 asym. stretch

---

Each simulation is initiated in the same way; vertical excitation of the ground-state vibrational wavefunction of a harmonic potential, i.e. a Gaussian of width in mass-frequency weighted coordinates, to the first excited singlet state (S₁). The individual GWPs were all centered at the FC point, but are distributed in momentum space. To represent an exact vertical excitation, only the wavepacket with no initial momentum is initially populated. The basis functions are 6-dimensional (3N − 6 normal modes), and a total of 49 GWPs was used for the two reference calculations. This number was chosen based on the formula N_GWP = n·2d + 1, where d is the dimensionality of the system and n is any integer. By using multiples of twice the number of DOFs, there should in principle be enough basis functions to describe motion along all normal modes in both directions. An additional GWP is then added to cover the FC region. This rationale, or a similar reasoning, has been applied in previous DD-vMCG studies,^52–54 as well as in quantum Ehrenfest calculations which use the same Gaussian basis functions.⁵⁵

Notice however, that this rule of thumb says nothing about whether or not the basis set is large enough for convergence. With respect to assessing the level of convergence of a DD-vMCG calculation, there are no rigorous criteria and it is therefore not our purpose to achieve fully converged results. That being said, based on preliminary convergence studies undertaken in Section S4 of the SI, the basis set used for our simulations is believed to be flexible enough to capture the important dynamics of the system. Since the molecule is not expected to undergo any large-scale geometric distortions (e.g. dissociation), the dynamics are limited to electronic phenomena mediated by normal mode oscillations. In the SI, we provide data showing the dependence of the population transfer with respect to the size of the nuclear basis, which indicates that no qualitative error is expected due to a lack of convergence. Attempts are made, referring to energy conservation, to argue for somewhat stronger quantitative assessments of the same. Unless otherwise stated, the results presented in the following section were obtained with a basis set of 49 GWPs – in line with our convergence analysis.

4 Results and discussion

This section is organised by first discussing the results of the 2 reference simulations. These are chosen based on the level of electronic structure theory used, and the stability of the propagation (in terms of total energy conservation for example, as shown in the SI). Calculations based on MS-CASPT2(12,10) and MRCI(10,10) potentials were selected (see Fig. 1 for electronic state energies at the FC point). The character of each state in the FC region is also given in Fig. 1. In terms of nomenclature, we shall use S₀, S₁, T₁, and T₂ labels when referring to specific states, but these will always be taken to correspond to the same electronic character as identified in the FC region. Since the PESs do not exhibit any major or important crossing points, this convention is appropriate and improves the clarity of the discussion to follow.

Fig. 1 Vertical excitation energies of TF (in eV) calculated with MS-CASPT2(12,10) and MRCI(10,10). Dashed lines are experimental values.⁵⁶

In the next sections, the resulting nuclear motion in both simulations is first described and a relevant cut of the PES is presented. Following this analysis, the diabatic state populations after excitation to S₁ are extracted for each of the four states, and the various possible population transfer processes are discussed, commenting on both IC, facilitated by NAC between S₀–S₁, and ISC, made possible by SOC elements between S₀–T₁ and S₁–T₂ with a magnitude around 160 cm⁻¹ (= 0.020 eV) for both pairs of states. We note that SOCs between S₀–T₂ and S₁–T₁ are negligible, and the corresponding transitions remain El-Sayed forbidden since no change in orbital angular momentum occurs.

Several additional DD-vMCG simulations were also carried out, based on different, but similar, quantum chemistry approaches. These are discussed in Section 4.3, where the sensitivity of direct dynamics simulations to the underlying potential is explored. Finally, a critical analysis of the observed differences between two non-adiabatic dynamics methods (TSH and DD-vMCG) is given in Section 4.4.

4.1 Nuclear motion

The most important nuclear DOFs during 1 ps of dynamics were identified based partially on the expectation values of the variance of each normal mode, as well as on previously observed trends. In addition to the C–S stretching mode (q₃), identified in prior studies, both reference calculations highlight the out-of-plane dihedral angle (q₁) as the two most active upon photoexcitation. The expectation values of these modes are plotted in Fig. 2 and 3, respectively, calculated from the full wavefunction (as given in eqn (1)) by projecting the expectation value of the position operator (〈Ψ|q̂|Ψ〉) onto each mode individually. The width of the overall wavepacket is then given by a Gaussian function where the variance equals that of the normal mode coordinate. In terms of units, we shall present all our results expressed as mass-frequency weighted normal modes, as were used during the simulations.

Fig. 2 Expectation value of the C–S bond stretching mode (q₃, mass-frequency weighted) during MS-CASPT2(12,10) (top) and MRCI(10,10) (bottom) dynamics. Signal centered on the expectation value of the position of the wavepacket, width given by its variance.

Fig. 3 Expectation value of the out-of-plane dihedral (q₁, mass-frequency weighted) during MS-CASPT2(12,10) (top) and MRCI(10,10) (bottom) dynamics. Signal centered on the expectation value of the position of the wavepacket, width given by its variance.

To see how nuclear motion affects the electronic states, a cut of the PES is presented in Fig. 4. The chosen coordinate is a linear combination of the q₁ and q₃ modes each having equal weights, thus corresponding to a 45° vector spanning the space in between the two orthogonal modes. Hence, the geometries along this coordinate are obtained by varying those two modes simultaneously, but keeping all others at q = 0.

Fig. 4 Top: PES cut along modes q₁ and q₃, generated from DD-vMCG dynamics at the MS-CASPT2(12,10) (solid) and MRCI(10,10) (dashed) levels of theory. Bottom: Energy difference between selected pairs of states. Units are given as mass-frequency weighted normal modes, specific values of the C–S bond length and dihedral angle are given above the figure.

4.2 Diabatic state populations

The diabatic state populations resulting from S₁ excitation are now analysed. These are easily extracted from state-specific components of the full vMCG wavefunction, whose magnitude yield the density P^diab_i on the ith diabatic state: P^diab_i(t) = |Ψ_i(t)|². The changes in population of all 4 states are plotted in Fig. 5, for the two reference simulations. The similarities and differences are described and explained below.

Fig. 5 Change in diabatic state population density during DD-vMCG simulations using MS-CASPT2(12,10) (top) and MRCI(10,10) (bottom) PESs.

4.2.1 IC: S₁ (¹nπ*) → S₀ (g.s.). While existing studies typically focus on ISC in photo-excited TF, it has nonetheless been identified as a typical case where the two NRR pathways of IC and ISC could compete.¹² The present study reveals interesting features related to IC in TF. Both reference calculations predict a certain amount of IC, with CASPT2 yielding about 8% and MRCI 5% after 1000 fs. In order to determine the mechanism behind this process, the minimum energy conical intersection (MECI) between the S₀–S₁ pair of states was located for both calculations, the result of which is summarised in Table 2.

Table 2 Energy and q₁ and q₃ normal mode components of the S₀–S₁ MECI, given both in mass-frequency weighted normal mode coordinates and in dihedral angle angle for q₁ and C–S bond length for q₃

	MS-CASPT2(12,10)	MRCI(10,10)
Energy	3.20 eV	3.99 eV
q1 (dihedral)	3.45 (114.9°)	3.42 (140.0°)
q3 (RC–S)	8.75 (2.16 Å)	15.17 (2.55 Å)

---

Based on the relative activity of the normal modes, both simulations agree on the mechanism of IC, being facilitated by a combination of the q₁ (out-of-plane bending) and q₃ (C–S stretching) normal modes. It might appear contradictory, at first glance, to find higher amounts of IC for CASPT2 compared to MRCI, in light of the former's higher S₀–S₁ energy gap (as found in the 1-D cuts and vertical excitation energies at the respective FC points). However, the nuclear motion along the aforementioned modes, as visualised in Fig. 2 and 3, reveals greater activity for CASPT2, especially after 500 fs when the IC rates start to differ. Looking at the MECIs described in Table 2, it is highly unlikely that any nuclear density reaches this point during the MRCI simulations, while the CASPT2 MECI could play a minor role. When looking at the nuclear density on the S₁ state (Fig. 6), we do indeed find limited instances where the wavepackets approach the MECI during the CASPT2 simulations, but not during the MRCI one.

Fig. 6 Nuclear density along modes q₁ and q₃ on the S₁ state at selected time-steps for both MS-CASPT2(12,10) (left), and MRCI(10,10) (right) simulations. The position of the MECI is indicated by the red dot. Timestamps were chosen based on the maximum displacement in the direction of the MECI.

The MRCI simulation suggests that TF remains remarkably photostable up to 1000 fs (and possibly beyond). Movement along the normal modes q₁ and q₃ remains constrained and consistent, accompanied by limited IC due to weak coupling between S₁ and high-lying vibrational levels of S₀, in agreement with experimental data.^57,58 We observe that the population of S₀ shows a faint and gradual increase for MRCI, suggesting that it is mediated by residual non-adiabatic coupling far away from regions of near-degeneracy (see Fig. S5 of the SI for a distribution of the NAC matrix elements). In contrast, IC for CASPT2 is sharper, driven instead by the movement of some of the nuclear wavepackets towards a lower-lying, more accessible MECI (3.20 eV), as seen in Fig. 6. However, it remains difficult to provide a definite quantitative prediction as to which of the two methods correctly captures the IC dynamics of TF, minimal as their differences are. Interestingly, the CASPT2 simulation, contrary to MRCI, suggests an alternate channel for ISC as discussed in the next section.

4.2.2 ISC: S₀ (g.s.) → T₁ (³nπ*). Here we give a brief discussion of the ISC channel made possible by SOC between the ground state and the lowest triplet state (El-Sayed allowed). The population of the ³nπ* state remains low in both reference calculations, but the CASPT2 calculation does predict a significantly higher amount of total ISC compared to MRCI, yielding around 5% after 1000 fs for the former and less than 1% for the latter. Additionally, the time scale and rate of population transfer differs in the two calculations (sharp increase after 500 fs vs. gradual rise), as was also observed – and discussed – with respect to IC. This is not surprising since the two processes are obviously connected; the ground state must first be populated before any ISC to T₁ can occur. Again, we explain this observation by referring to the differences in nuclear motion along normal modes q₁ and q₃, and the energies of the MECI for both simulations (Table 2). The bottom panels of Fig. 4 reveal that the energy difference between S₀–T₁ does indeed decrease sharply along the combined coordinate. Even though MRCI predicts a smaller energy gap than CASPT2, the magnitude of the nuclear motion is found to be larger during simulations based on CASPT2 surfaces (Fig. 2 and 3), thus extending into regions of the configuration space where the energy gap is smaller than in the regions explored by MRCI dynamics. Importantly though, it is noted that in both simulations the level of T₁ population eventually exceeds that of T₂. This observation should be expected in any simulation which is extended to long times, say >500 fs, since the increasing S₀ population will continue to transfer some fraction of the density to T₁, while the rate of T₂ population remains negligible.

We also note that although ΔE between S₁–T₂ is initially smaller than S₀–T₁, it does not decrease as rapidly with respect to q₁ and/or q₃. This explains why no significant T₂ population is observed, even at later times, despite the SOC strength being essentially equal between the S₀–T₁ and S₁–T₂ pairs.

With respect to previous non-adiabatic dynamics studies,^11,14–16 it seems to be the consensus that no significant T₁ population is expected within 500 fs. We agree with this assertion in that specific time window, but if one accepts that the IC channel becomes important at longer time scales, the possibility of this alternative ISC pathway should not be ignored. Since those investigations do not discuss IC, it is unclear whether or not the possibility for this process, and by extension ISC to T₁, was included in the simulations or not.

4.2.3 ISC: S₁ (¹nπ*) → T₂ (³ππ*). Finally, we comment on population transfer to the highest lying triplet state. The reference simulations both agree on the fact that no significant density builds up in this state during 1 ps of dynamics. This observation also matches experimental evidence that no significant ISC is expected on sub-picosecond timescales.⁵⁶ There are however periodic features showing alternating population and de-population of the triplet, but always staying below 1%. In order to rationalise this observation, we start by pointing out that although SOC between ¹nπ* and ³ππ* is non-zero at the FC point, it is still relatively weak and the energy gap is too large (>1 eV) to allow any population transfer at this geometry. It follows then, that any ISC observed between the states must be a result of some nuclear motion modulating this singlet–triplet gap. The nature of this motion is inferred by looking at Fig. 2 and noticing that the coherent oscillations (persisting during the first 800 fs of CASPT2 dynamics, and during the full 1000 fs of the MRCI simulation) match the frequency of the small rises observed in the T₂ population. The period is found to be 44 fs in both cases, in good agreement with the experimentally measured value of 41 fs.⁵⁶ Because of this connection, we propose a mechanism where, as the C–S bond lengthens, a very small amount of density is transferred from the singlet to the triplet state, reaching a maximum when the bond length is at its longest (〈R_C–S〉 = 1.93 Å). Then, the density is transferred back to the singlet as the bond starts to shorten. The bottom part of Fig. 4 also confirms that the C–S bond stretch coordinate does indeed bring the S₁ and T₂ states closer together.

4.3 Sensitivity to electronic structure theory

Within the Born–Oppenheimer framework, the dynamics of a molecule are described by the movement of nuclear wavepackets over a set of potential energy surfaces (PESs), obtained from an appropriate level of electronic structure theory. It is expected therefore, that simulating non-adiabatic dynamics fundamentally depends on the choice of electronic structure method vis-à-vis the PESs. This sensitivity is further expressed for direct-dynamics methods, where the PESs are calculated on-the-fly, instead of being fit and regularised over a global grid. It is also expected that the dynamics of certain systems should be more sensitive than others, in which case the level of electronic structure must be carefully determined.

Previous studies with trajectory surface hopping (TSH)^11,14,15 highlight TF as one such case where the dynamics appear vastly different across different electronic structure methods, despite there being a general agreement among them at equilibrium, in terms of optimised geometries and excitation energies. However, the present set of DD-vMCG calculations, contrary to TSH, shows no major differences in terms of state populations up to 500 fs of simulation. Each of the three different active spaces employed [(12,10), (10,10), (10,6)], and the dynamically correlated MS-CASPT2 and MRCI variations predict that TF remains photostable up to 500 fs, showing very little to no IC or ISC, and no prominent long-range motion, with differences only apparent at longer time-scales of 1 ps. The varying sensitivities to the underlying PESs found in these studies could be attributed to the different dynamics methods employed, and/or the different state representations utilized to describe SOC therein. A more elaborate comparison follows in the subsequent section. Here, we first discuss the other DD-vMCG simulations and compare them to our references (results summarised in Table 3, energy conservation and population plots given in the SI).

Table 3 A summary of the dynamics calculations performed. State populations are indicated as (final, at 500 fs). MC: Molcas,⁵⁹ MP: Molpro⁶⁰

Method	PES		Dynamics		Population (%)
	Energy gaps	Freq.	NGWP	Time (fs)	S0	T1	T2
REFERENCE:
MS-CASPT2(12,10) (MC)	✓	✓	49	1000	8,1	5,0	≈0,0
MRCI(10,10) (MP)	✓	✓	49	1000	5,1.5	<1,0	≈0,0
MRCI(12,10) (MP)	✓	✓	37	800	8.5,<3	<1,0	≈0,0
CASSCF(12,10) (MP)	Low	Low	25	1000	4.5,1.5	<1,0	<1,0
CASSCF(10,10) (MP)	Low	Low	37	1000	4.5,<1	<1,0	≈0,0
CASSCF(10,6) (MC)	Low S1–T2	✓	49	1000	10,1	2,<1	2,<1
CASSCF(10,6)(-Opt) (MC)	Low S1–T2	✓	49	950	15,3	<5,<1	<5,1

---

As the smallest of the active spaces proposed, dynamics on the CASSCF(10,6) surfaces predict reduced fluorescence yields for TF at longer time-scales. Although ISC to the El-Sayed allowed T₂ state remains low, unusually high rates of IC to the ground state are observed. A second set of calculations, using an optimised geometry and frequencies similar to those of previous TSH studies, predicts even higher rates of IC and ISC. Comparing these to our reference calculations, we conclude that the dynamics of TF are well described with the CASSCF(10,6) surfaces up to 500 fs, despite the low S₁–T₂ energy gap. However, the smaller active space does not extend to regions of the nuclear configuration space further displaced from equilibrium geometry – one suspects the exclusion of the C–H σ* orbital from the active space affects the electronic structure at stretched nuclear geometries. Consequently, the calculations become replete with artifacts at longer time-scales of 1 ps as the nuclear wavepackets spread out. In particular, normal modes q₅ and q₆, associated with C–H stretching, show tendencies of long-range motion. Coupled with the generally low S₁–T₂ energy gap, oscillations in T₂ population (beyond 500 fs) give way to a more pronounced, monotonic increase (see Fig. S6 in the SI).

Dynamics over the dynamically-correlated MRCI PESs, based on the (10,10) active space, forms one of the two references in the present study of photo-excited TF. Although the active space remains well-behaved for the entirety of the simulations, the CASSCF(10,10) dynamics reveal a number of artifacts. Unlike the slow and gradual IC of the MRCI reference, the CASSCF(10,10) calculations show very little IC up to 850 fs, followed by a sudden and erratic transfer of population to the S₀ state thereafter, uncharacteristic for a generally fluorescent molecule. This is possibly due to the accumulation of numerical errors at long simulation time-scales, as reflected in the loss of total energy conservation. While long-time propagation in quantum dynamics simulations is bound to introduce such numerical artifacts, these are found not to be quite so prominent in this case – energy is reasonably conserved over 1 ps (see Fig. S7 in the SI). Instead, the CASSCF(10,10) surfaces, obtained from Molpro, systematically underestimate the S₁–S₀ energy gap (inflating IC), as well as the optimised frequencies of the IC-mediating q₁ mode (deflating IC), the latter utilized in the definition of the kinetic energy operator for the nuclear wavepackets. Thus, despite there not being major differences in the final state populations with the MRCI reference, the CASSCF(10,10) PESs do not generate correct dynamics. This is resolved upon introducing dynamic correlation in the description of electronic structure with MRCI, and the resulting dynamics appear well-behaved with respect to energy conservation metrics, even at longer times, and in good agreement with experimental data. Compared to the more intuitive (12,10) variant, it is interesting to note that the (10,10) active space with MRCI provides well-behaved PESs.

As the most extensive of the lot, the (12,10) active space remains stable over the entire range of nuclear geometries explored, with the MS-CASPT2 calculations serving as a reference. Although the CASPT2 (from Molcas) and MRCI PESs (from Molpro) show the right features, the CASSCF(12,10) PESs (from Molpro) underestimate energy gaps and optimised frequencies of key normal modes – quite like CASSCF(10,10), albeit to a lesser extent. This affects the resulting dynamics as the state populations lack the generally smooth and gradual features found in the reference calculations. Similarly, the MRCI dynamics predict a premature loss (around 500 fs) of the periodic oscillations of the T₂ state associated with the C–S stretch mode (q₃) followed by unusual behaviour along the C–H stretch (symmetric and anti-symmetric) coordinates around 800 fs, when the calculations are terminated. For both CAS and MRCI calculations, these anomalies arise due to the dissociative behaviour of some spurious wavepackets along q₅ and q₆ modes (C–H stretch) – although the molecule, on average, does not seem to develop any long-range motion. Normal modes q₅ and q₆, skewed by these spurious wavepackets, appear more active in these calculations compared to the two references – this is possibly due to sub-optimal coupling between GWPs in the nuclear basis, which could only be represented with 25 GWPs (for CAS) and 37 GWPs (for MRCI; only up to 800 fs), due to numerical complications associated with DD-vMCG (discussed next). This is despite the general stability of the (12,10) active space, as reflected in the MS-CASPT2 results. However, the aforementioned long-time propagation errors are somewhat apparent for the MS-CASPT2 simulations (see Fig. S7 in the SI) – for an on-the-fly calculation, the loss of energy conservation is exacerbated by the uneven nature of the PESs (since these are not smooth, analytic functions) under LHA, alongside errors associated with GWP propagation itself. While one could, in principle, improve the quality of the PESs by calculating more ab initio points, this is not expected to majorly affect the outcome of the simulation – in terms of final state populations, for instance – if already within an acceptable convergence threshold. The SI contains a more detailed convergence analysis. Results from the different dynamics simulations are tabulated in Table 3.

We conclude this section by discussing one of the well-known complications associated with the DD-vMCG method, namely, the linear dependence of the Gaussian wavepackets.^61–64 This issue results in a potential practical roadblock when performing DD-vMCG simulations. While it is generally desirable to represent the nuclear degrees of freedom with a sufficiently large number of Gaussian wavepackets – this is to ensure variational convergence to the true result – it is often the case, due to the nature of the underlying PESs, that increasing the number of wavepackets beyond a certain threshold produces artifacts due to the associated linear dependencies. The wavepackets are found to strongly overlap, typically in bound regions of the nuclear configuration space, where they cannot diffuse freely due to the presence of barriers. For the present study, while the two references and the CASSCF(10,6) calculations could be extended to 49 GWPs (in accordance with the heuristic N_GWP = n·2d + 1), the other calculations begin to show irregularities when increasing the number of wavepackets. Therefore, the SCF (10,10) and MRCI (12,10) dynamics could only be performed with 37 GWPs, with an even lower 25 GWPs for SCF (12,10), resulting in dissociative artifacts for the latter two. In the event the calculations could be extended with more wavepackets, the state populations are expected to change and move closer to convergence. This is illustrated in the convergence plots for the MS-CASPT2 calculations (see the SI).

4.4 Sensitivity to QD method

Here we make a direct comparison between a DD-vMCG calculation and a TSH simulation taken from ref. 15. The reason for this comparison is not to comment on the ultimate accuracy (in fact, both simulations are in disagreement with experiment), but rather to examine some of the similarities and differences between the two methods. In accordance with the aims laid out in the QD method benchmarking roadmap,²⁷ this is believed to be a useful exercise. We also point out that previous studies exist which specifically compare TSH and (DD)-vMCG, but these focus on systems exhibiting rapid and pronounced population transfer.^65,66 The slower and more subtle effects in TF thus offer a new dimension to the comparison of these two methods.

One important difficulty related to comparing dynamical methods lies in the fact that one must ensure that the underlying potential is identical, so as to only evaluate differences in nuclear propagation. We have thus chosen to follow the same procedure as outlined in ref. 15: geometry optimisation using SA-CASSCF(10,6) with the cc-pVDZ basis set, and subsequently using that method for on-the-fly dynamics. Following the propagation, the state populations were calculated and are summarised in Table 4 for both the CASSCF(10,6) simulations performed with DD-vMCG in this work, and those performed with TSH (implemented in SHARC^67,68) in ref. 15.

Table 4 Comparison of DD-vMCG and TSH results after 500 fs, using SA-CASSCF(10,6)

	DD-vMCG (this work)	TSH^15
S0	3%	—
T1	<1%	≈0%
T2	1%	5%

---

Starting with IC back to the ground state (S₀), no comparison can be made since the ground state population is not reported for the TSH result. It is thus inferred that it was not found to play an important role during the dynamics, and as a result, the population of T₁ is also not expected. Therefore we focus on the T₂ population, which is small in both simulations but not zero (as expected). The qualitative agreement between the oscillations observed in both simulations seems to indicate that the coupling between S₁–T₂ has a very similar effect, despite the use of different spin representations. This strengthens our assertion that using a spin-diabatic formalism can be appropriate in this case. Quantitatively speaking, the relative magnitude of the overestimation does differ between the two QD methods, indicating that the different formalisms used to propagate the nuclei are also playing a role here. Although NRR remains marginal in either simulation, we may detect hints of a competition between IC and ISC channels – one could argue it is IC from S₁ to S₀, in turn followed by ISC from S₀ to T₁, which limits the extent of ISC from S₁ to T₂ for the DD-vMCG simulation. We wish to evaluate this aspect as we believe it is of general interest to the QD community to better understand how exactly various methods differ from one another. In what follows, a selection is made of several possible factors that may affect the final populations of the states in the two calculations.

One of the most obvious differences between DD-vMCG and TSH is the coupled nature of the GWPs in the former, versus the independent trajectories of the latter. A result of this difference in nuclear basis set is that even though the same electronic structure method is used, the configuration space explored by the two approaches will not be identical, and therefore neither are the on-the-fly PESs. Since the photophysical processes in TF are somewhat sensitive to the underlying potential, the difference in explored geometries may be entirely responsible for the discrepancies in Table 4. One way to check this hypothesis would be to run TSH dynamics on the database generated from the DD-vMCG simulation, to see if the resulting dynamics and state populations bear more of a resemblance to the presented DD-vMCG results. Although not currently implemented for both singlet and triplet states, this capability is expected to be added to Quantics in the near future.

Instead, we investigate the importance of coupling between basis functions during the propagation. Using the database generated by CASSCF(10,6) DD-vMCG, a second simulation was run where the variational coupling between GWPs was turned off (so-called classical MCG, clMCG^23,69). This had virtually no effect on the final state populations, indicating that classical trajectories and coupled wavepackets may give very similar results for TF when run on the exact same PES. To test this hypothesis more strongly, we also rule out the effects of convergence, or rather the lack thereof, by investigating the effect of basis size. Four additional DD-vMCG and DD-clMCG simulations were performed with n = 1, 2, 3, and 4 in the formula N_GWP = n·2d + 1, i.e. 13, 25, 37, and 49 GWPs respectively. In all cases, the T₂ populations were not significantly affected.

Next, the generation of initial conditions, i.e. the choice of Ψ(t = 0), is considered. Again, TSH and vMCG differ in this respect but they do also both rely on similar approximations, namely the assumption that (1) an infinitely short excitation pulse is used, and (2) the wavepacket is constructed from the system in its electronic and vibrational ground state (a valid assumption at 0 K). The difference between the methods is that in vMCG, the vibrational wavefunction of the ground state (i.e. a Gaussian, in the harmonic approximation) is projected onto an excited state, whereas in TSH, the initial momenta and positions for the independent trajectories are sampled from a distribution constructed from the ground state geometry and its normal modes (e.g. using a harmonic Wigner distribution⁷⁰). Clearly, direct comparison of two methods based on different initial conditions is not straightforward, as was also identified in the benchmarking roadmap. Furthermore, achieving consistent initial conditions between two methods based on fundamentally different formalisms (trajectories vs. wavefunctions, in this case) remains an open problem. In an effort to address this question, a DD-vMCG simulation was initialised where each of the 49 GWPs were assigned initial momenta and positions according to a Wigner distribution (individual components of the distribution are shown in the SI). However, this calculation suffered from similar numerical issues to those discussed in the previous section, only starting much earlier at about 200 fs, indicating that these conditions are not suitable for vMCG. We do note that no population transfer had occurred in that time.

Related to the discussion of initialisation, the question of state representation should be mentioned when directly comparing methods. As stated previously, our simulations employ the diabatic representation, whereas SHARC uses a “diagonal” representation of the Hamiltonian^71,72 which is more closely related to the adiabatic picture, in the sense that the potential is constructed from the eigenvalues of the total electronic Hamiltonian. As a result, the PESs of the TSH simulations are not necessarily spin-pure, and triplet state population can be observed without any hops directly to a (dominantly) triplet state. In the present case however, the build up of about 5% T₂ population is too high to be attributed to this effect and we thus conclude that the rather subtle difference in state representation is not likely to cause the discrepancy between the two results.

With regards to the observation of IC to S₀, there are striking differences in the way NAC matrix elements (NACMEs) are described and transformed in the two propagation algorithms – while DD-vMCG makes use of so-called propagation diabatisation to transform the ab initio NACMEs to a diabatic frame, the particular implementation of TSH in this comparison uses local diabatisation⁷³ to extract hopping probabilities across adiabatic electronic states, without explicit computation of NACMEs. Since IC in TF is not so much mediated by nuclear motion in the vicinity of conical intersections, but rather by residual couplings spread variously across the nuclear configuration space (see the SI for distribution of NACMEs), the sensitivity of the chosen diabatisation scheme to these couplings could also have an influence on the dynamics. This, however, is non-trivial to quantify.

Based on the above discussion, we hypothesise that the most likely factor causing the slight discrepancies between DD-vMCG and TSH is related to how the phase-space of the molecule is explored during on-the-fly dynamics. This indicates that both the differences in initial conditions, and the coupled versus independent nature of the nuclear basis may play an important role in this exploration. It is however difficult to provide conclusive evidence at this time, and furthermore, this conclusion is not easily extendable to the general case. To better quantify these differences, it may be useful to define concrete measures for how different regions of the PES are explored during dynamics. At the moment, a proper framework to discuss these effects is lacking, but it is expected that community-driven efforts towards suitable benchmarks for QD methods (e.g. the ideas included in the aforementioned roadmap) will help put these results in perspective and allow for a more thorough discussion.

5 Conclusion

Following the successful implementation of spin–orbit coupling in fully coupled direct wavepacket dynamics, the present study forms the first instance of the DD-vMCG method, in the Quantics suite of programs, utilized to describe photo-excited dynamics over a singlet–triplet manifold. Specifically, we quantify the extent of non-radiative intersystem crossing to different triplet states of TF, following photo-excitation to the S₁ state. Although TF is found to remain generally fluorescent – with no prominent NRR channels, and photostable – with no significant long-range motion, up to the picosecond time-scale, in agreement with existing theoretical and experimental studies, the dynamics reveal an interesting competition between the two nonradiative pathways of IC and ISC. Our two reference calculations – on the MRCI(10,10) and MS-CASPT2(12,10) PESs – predict a small but consistent presence of IC(5–8%), and the accompanying possibility of ISC to T₁, mediated by a combination of the C–S stretching (q₃) and out-of-plane bending (q₁) modes. The photophysics of TF is found to be sensitive to the level of electronic structure used to derive the underlying PESs, based on differences in state populations at longer time-scales of simulation. However, the nature and extent of this sensitivity could be subject to the dynamics method chosen. In line with the ideas laid out in the roadmap for QD benchmarks, we compared our findings with existing TSH studies that utilize the same level of electronic structure theory. A preliminary attempt is made to explain the differences, and some of the difficulties associated with benchmarking are described. As these difficulties are overcome, the subtle features of the excited state dynamics of TF may offer fresh prospects for future benchmark studies.

Author contributions

Swagato Saha: methodology; investigation; formal analysis; visualisation; writing – original draft; writing – review & editing. Léon Cigrang: methodology; investigation; formal analysis; visualisation; writing – original draft; writing – review & editing. Graham Worth: writing – review & editing (equal); project administration; supervision; funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information (SI): including the active orbitals, coordinates, demonstration of convergence and initial conditions. See DOI: https://doi.org/10.1039/d5cp04004c.

The Quantics code used for the calculations is open source and available on request to the authors (https://www.chem.ucl.ac.uk/quantics). The input files for the calculations and the databases produced by the DD-vMCG calculations are available from the UCL data repository (https://rdr.ucl.ac.uk) at the DOI: https://doi.org/10.5522/04/30366454.

Acknowledgements

This work was performed as part of the EPSRC programme grant COSMOS EP/X026973/1. SS also acknowledges a PhD Studentship from the UCL EPSRC Doctoral Training Account EP/W524335/1.

Notes and references

1. C. M. Marian, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 187–203.
2. J. Li and S. A. Lopez, Computational Chemistry for Photochemical Reactions, Elsevier, 2024, pp. 658–698.
3. A. Gilbert and J. Baggott, Essentials of Molecular Photochemistry, Blackwell Scientific Publications, 1991.
4. G. N. Lewis and M. Kasha, J. Am. Chem. Soc., 1945, 67, 994–1003.
5. M. Kasha, Chem. Rev., 1947, 41, 401–419.
6. J. Zhao, W. Wu, J. Sun and S. Guo, Chem. Soc. Rev., 2013, 42, 5323.
7. K. Miyata, F. S. Conrad-Burton, F. L. Geyer and X.-Y. Zhu, Chem. Rev., 2019, 119, 4261–4292.
8. X.-K. Chen, D. Kim and J.-L. Brédas, Acc. Chem. Res., 2018, 51, 2215–2224.
9. R. S. Minns, D. S. N. Parker, T. J. Penfold, G. A. Worth and H. H. Fielding, Phys. Chem. Chem. Phys., 2010, 12, 15607–15615.
10. C. Reichardt, R. A. Vogt and C. E. Crespo-Hernández, J. Chem. Phys., 2009, 131, 224518.
11. L. Zhang, Y. Shu, S. Sun and D. G. Truhlar, J. Chem. Phys., 2021, 154, 094310.
12. Y. Guan, C. Xie, H. Guo and D. R. Yarkony, J. Chem. Theory Comput., 2023, 19, 6414–6424.
13. D. C. Moule and E. C. Lim, J. Phys. Chem. A, 2002, 106, 3072–3076.
14. P. B. Calio, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2022, 18, 614–622.
15. S. Mai, A. J. Atkins, F. Plasser and L. González, J. Chem. Theory Comput., 2019, 15, 3470–3480.
16. B. F. Curchod, C. Rauer, P. Marquetand, L. González and T. J. Martínez, J. Chem. Phys., 2016, 144, 101102–101106.
17. G. W. Richings and S. Habershon, J. Phys. Chem. A, 2020, 124, 9299–9313.
18. H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.
19. U. Manthe, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199–3213.
20. M. H. Beck, A. Jäckle, G. A. Worth and H.-D. Meyer, Phys. Rep., 2000, 324, 1–105.
21. H.-D. Meyer and G. A. Worth, Theor. Chem. Acc., 2003, 109, 251–267.
22. G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt and B. Lasorne, Int. Rev. Phys. Chem., 2015, 34, 269–308.
23. K. E. Spinlove, M. Vacher, M. Bearpark, M. A. Robb and G. A. Worth, Chem. Phys., 2017, 482, 52–63.
24. G. A. Worth, M. A. Robb and I. Burghardt, Faraday Discuss., 2004, 127, 307–323.
25. G. W. Richings and G. A. Worth, Chem. Phys. Lett., 2017, 683, 606–612.
26. G. A. Worth, Comput. Phys. Commun., 2020, 248, 107040.
27. L. 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ć, E. C. Eklund, S. Fernandez-Alberti, A. Freibert, L. González, G. Granucci, F. J. Hernández, J. Hernández-Rodríguez, A. Jain, J. Janoš, I. Kassal, A. Kirrander, Z. Lan, H. R. Larsson, D. Lauvergnat, B. L. Dé, Y. Lee, N. T. Maitra, S. K. Min, D. Peláez, D. Picconi, U. Raucci, Z. Qiu, P. Robertson, E. S. Gil, M. Sapunar, P. Schürger, P. Sinnott, S. Tretiak, A. Tikku, P. Vindel-Zandbergen, G. A. Worth, F. Agostini, S. Gómez, L. M. Ibele and A. Prlj, J. Phys. Chem. A, 2025, 129, 7023–7050.
28. G. A. Worth and L. S. Cederbaum, Annu. Rev. Phys. Chem., 2004, 55, 127–158.
29. D. R. Yarkony, Rev. Mod. Phys., 1996, 68, 985–1013.
30. A. Migani and M. Olivucci, Conical Intersections and Organic Reaction Mechanisms, World Scientific, 2004, ch. 6, pp. 271–320.
31. G. A. Meek and B. G. Levine, J. Chem. Phys., 2016, 144, 184109.
32. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat and F. Agostini, Phys. Rev. Lett., 2020, 124, 033001.
33. G. Granucci, M. Persico and G. Spighi, J. Chem. Phys., 2012, 137, 22A501.
34. T. J. Frankcombe, M. A. Collins and G. A. Worth, Chem. Phys. Lett., 2010, 489, 242–247.
35. T. J. Frankcombe, J. Chem. Phys., 2014, 140, 114106–114108.
36. G. W. Richings and G. A. Worth, J. Phys. Chem. A, 2015, 119, 12457–12470.
37. M. Baer, Chem. Phys. Lett., 1975, 35, 112–118.
38. C. Sanz-Sanz and G. A. Worth, Phys. Chem. Chem. Phys., 2019, 21, 14429–14439.
39. F. Talotta, S. Morisset, N. Rougeau, D. Lauvergnat and F. Agostini, J. Chem. Theory Comput., 2020, 16, 4833–4848.
40. B. O. Roos, Adv. Chem. Phys., 1987, 87, 399–466.
41. H.-J. Werner and W. Meyer, J. Chem. Phys., 1981, 74, 5794–5801.
42. K. Anderson, P.-Å. Malmqvist, B. O. Roos, A. J. Sadlej and K. Wolinski, J. Phys. Chem., 1990, 94, 5483–5488.
43. J. Finley, P.-A. Malmqvist, B. O. Roos and L. Serrano-andrés, Chem. Phys. Lett., 1998, 288, 299–306.
44. F. Plasser and H. Lischka, Multi-Reference Configuration Interaction, Wiley, 2020, ch. 9, pp. 277–297.
45. J. Janoš and P. Slavíček, J. Chem. Theory Comput., 2023, 19, 8273–8284.
46. D. Bellshaw, R. S. Minns and A. Kirrander, Phys. Chem. Chem. Phys., 2019, 21, 14226–14237.
47. D. Tuna, D. Lefrancois, L. Wolański, S. Gozem, I. Schapiro, T. Andruniów, A. Dreuw and M. Olivucci, J. Chem. Theory Comput., 2015, 11, 5758–5781.
48. T. V. Papineau, D. Jacquemin and M. Vacher, J. Phys. Chem. Lett., 2024, 15, 636–643.
49. L. M. Ibele, A. Memhood, B. G. Levine and D. Avagliano, J. Chem. Theory Comput., 2024, 20, 8140–8151.
50. S. Mukherjee, R. S. Mattos, J. M. Toldo, H. Lischka and M. Barbatti, J. Chem. Phys., 2024, 160, 154306.
51. G. A. Worth, K. Giri, G. W. Richings, I. Burghardt, M. H. Beck, A. Jäckle and H.-D. Meyer, Quantics package, Version 2.2, 2024.
52. L. L. E. Cigrang and G. A. Worth, J. Phys. Chem. A, 2024, 128, 7546–7557.
53. O. Bennett, A. Freibert, K. E. Spinlove and G. A. Worth, J. Chem. Phys., 2024, 160, 174305.
54. K. E. Spinlove, G. W. Richings, M. A. Robb and G. A. Worth, Faraday Discuss., 2018, 212, 191–215.
55. T. Tran, G. A. Worth and M. A. Robb, J. Comput. Chem., 2025, 46, e70028.
56. D. J. Clouthier and D. A. Ramsay, Annu. Rev. Phys. Chem., 1983, 34, 31–58.
57. J. R. Dunlop and D. J. Clouthier, J. Chem. Phys., 1990, 90, 6371–6386.
58. D. J. Clouthier, G. Huang and A. J. Merer, J. Chem. Phys., 1992, 97, 1630–1637.
59. G. Li Manni, I. F. Galván, A. Alavi, F. Aleotti, F. Aquilante, J. Autschbach, D. Avagliano, A. Baiardi, J. J. Bao, S. Battaglia, L. Birnoschi, A. Blanco-González, S. I. Bokarev, R. Broer, R. Cacciari, P. B. Calio, R. K. Carlson, R. Carvalho Couto, L. Cerdán, L. F. Chibotaru, N. F. Chilton, J. R. Church, I. Conti, S. Coriani, J. Cuéllar-Zuquin, R. E. Daoud, N. Dattani, P. Decleva, C. de Graaf, M. G. Delcey, L. De Vico, W. Dobrautz, S. S. Dong, R. Feng, N. Ferré, M. Filatov(Gulak), L. Gagliardi, M. Garavelli, L. González, Y. Guan, M. Guo, M. R. Hennefarth, M. R. Hermes, C. E. Hoyer, M. Huix-Rotllant, V. K. Jaiswal, A. Kaiser, D. S. Kaliakin, M. Khamesian, D. S. King, V. Kochetov, M. Krośnicki, A. A. Kumaar, E. D. Larsson, S. Lehtola, M.-B. Lepetit, H. Lischka, P. López Ríos, M. Lundberg, D. Ma, S. Mai, P. Marquetand, I. C. D. Merritt, F. Montorsi, M. Märchen, A. Nenov, V. H. A. Nguyen, Y. Nishimoto, M. S. Oakley, M. Olivucci, M. Oppel, D. Padula, R. Pandharkar, Q. M. Phung, F. Plasser, G. Raggi, E. Rebolini, M. Reiher, I. Rivalta, D. Roca-Sanjuán, T. Romig, A. A. Safari, A. Sánchez-Mansilla, A. M. Sand, I. Schapiro, T. R. Scott, J. Segarra-Martí, F. Segatta, D.-C. Sergentu, P. Sharma, R. Shepard, Y. Shu, J. K. Staab, T. P. Straatsma, L. K. Sørensen, B. N. C. Tenorio, D. G. Truhlar, L. Ungur, M. Vacher, V. Veryazov, T. A. Voß, O. Weser, D. Wu, X. Yang, D. Yarkony, C. Zhou, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2023, 19, 6933–6991.
60. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O'Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson and M. Wang, MOLPRO, version 2015.1, a package of ab initio programs, 2015, https://www.molpro.net.
61. S. Sawada, R. Heather, B. Jackson and H. Methiu, J. Chem. Phys., 1985, 83, 3009–3027.
62. S. Sawada and H. Methiu, J. Chem. Phys., 1986, 84, 227–238.
63. I. Burghardt, K. Giri and G. A. Worth, J. Chem. Phys., 2008, 129, 174104–174114.
64. S. Habershon, J. Chem. Phys., 2012, 136, 14108–14109.
65. J. Coonjobeeharry, K. E. Spinlove, C. Sanz Sanz, M. Sapunar, N. Došlić and G. A. Worth, Philos. Trans. R. Soc., A, 2022, 380, 20200386.
66. S. Gómez, E. Spinlove and G. Worth, Phys. Chem. Chem. Phys., 2024, 26, 1829–1844.
67. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola and L. González, J. Chem. Theory Comput., 2011, 7, 1253–1258.
68. S. Mai, P. Marquetand and L. González, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2018, 8, e1370.
69. L. L. E. Cigrang, J. A. Green, S. Gómez, J. Cerezo, R. Improta, G. Prampolini, F. Santoro and G. A. Worth, J. Chem. Phys., 2024, 160, 174120.
70. J. P. Dahl and M. Springborg, J. Chem. Phys., 1988, 88, 4535–4547.
71. S. Mai, P. Marquetand and L. González, J. Chem. Phys., 2014, 140, 204302.
72. F. Plasser, S. Gómez, M. F. S. J. Menger, S. Mai and L. González, Phys. Chem. Chem. Phys., 2019, 21, 57–69.
73. G. Granucci, M. Persico and A. Toniolo, J. Chem. Phys., 2001, 114, 10608.

---