# Ultrafast dynamics of photo-excited 2-thiopyridone: Theoretical insights into triplet state population and proton transfer pathways.

Jesper Norell, Michael Odelius, Morgane Vacher

► To cite this version:

Jesper Norell, Michael Odelius, Morgane Vacher. Ultrafast dynamics of photo-excited 2-thiopyridone: Theoretical insights into triplet state population and proton transfer pathways. Structural Dynamics, 2020, 7 (2), pp.024101. 10.1063/1.5143228. hal-03018929

HAL Id: hal-03018929

https://hal.science/hal-03018929

Submitted on 3 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ultrafast dynamics of photo-excited 2-thiopyridone: Theoretical insights into triplet state population and proton transfer pathways

EP

Cite as: Struct. Dyn. **7**, 024101 (2020); doi: 10.1063/1.5143228

Submitted: 20 December 2019 · Accepted: 29 February 2020 ·

Published Online: 17 March 2020

<figure><img src="image_1.png" /></figure>

View Online

<figure><img src="image_2.png" /></figure>

Export Citation

<figure><img src="image_3.png" /></figure>

CrossMark

Jesper Norell,<fnref n="1" />,ᵃ⁾ Michael Odelius,<fnref n="1" />,ᵇ⁾ and Morgane Vacher<fnref n="2" />,<fnref n="3" />,ᶜ⁾

## AFFILIATIONS

## ABSTRACT

Ultrafast non-adiabatic dynamics of the small heteroaromatic compound 2-thiopyridone has been studied with surface hopping simulations based on multi-configurational quantum chemistry. Initial excitation of the bright $S_2(\pi, \pi^*)$ state is found to promptly relax to $S_1(n, \pi^*)$ through in-plane motion. The subsequent dynamics are oppositely driven by out-of-plane motion, which results in both complex population transfers among all the available states and intersystem crossing predominantly through the “El-Sayed forbidden” $S_1(n, \pi^*)$ to $T_2(n, \pi^*)$ channel, through significant mixing of electronic excitation characters. Despite this complexity, the femto- to picosecond triplet population, expected from several spectroscopic measurements, is well described as a simple exponential decay of the singlet state manifold. No proton transfer is found in the reported trajectories, but two mechanisms for its possible mediation in previously reported experiments are proposed based on the observed structural dynamics: (i) ultrafast intra-molecular transfer driven by the initially coherent in-plane motion and (ii) inter-molecular solvent-mediated transfer driven by the out-of-plane modes that dominate the later motion.

© 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5143228

## I. INTRODUCTION

Due to the importance of light induced processes in biomolecules, in general, and in nucleic acids, in particular, the photochemistry of heteroaromatic compounds has been the subject of extensive investigations.<fnref n="1" /> Among these, thiosubstituted (organosulfuric) compounds have gathered increasing interest from ultrafast science due to their wide range of applications and rich excited-state dynamics.<fnref n="2" />⁻¹⁷ One such compound, 2-thiopyridone (2-TP, see Fig. 2), has primarily been studied as a compact model system for excited-state proton transfer (ESPT),¹⁸⁻²⁴ which is related to, for instance, photo-protection mechanisms of DNA²⁵ and pigmentation of the human skin.²⁶

While it is long established that the nitrogen protonated 2-TP can tautomerize into sulfur protonated 2-mercaptopopyridine upon a change in the chemical environment (gas phase, solid, or solvated),²⁷⁻³¹ it is only recently that time-resolved measurements¹⁸,²⁰⁻²² have been

used to explicitly investigate the possibility for the corresponding UV-pumped process. In the measurement with the highest temporal resolution, traces of N-H dissociation were found below picosecond (ps) timescales using N1s resonant inelastic x-ray scattering (RIXS), but, due to limited statistics, the electronic states involved in the process could not be resolved.²⁰ Two pico- to nanosecond (ns) measurements, performed using N1s²² and S1s²¹ near-edge x-ray absorption (NEXAFS), however, instead find the ESPT to be a nanosecond process, consistent with previous vibrational measurements.¹⁸ In particular, our recent study,²² combining time-resolved N1s NEXAFS with spectrum simulations of valence-excited states, provided strong evidence for a dominant reaction pathway on these timescales, shared for $S_2$ and $S_4$ initial excitation, where ESPT results from nanosecond decay of the $T_1$ state. Notably, none of the studies have been able to address whether the ESPT occurs intra-molecularly or inter-molecularly through solvent-mediated formation of a deprotonated anion $2-TP⁻$.

<footnote n="1">¹Department of Physics, AlbaNova University Center, Stockholm University, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="2">²Department of Chemistry-Ångström Laboratory, Uppsala University, 75121 Uppsala, Sweden</footnote>

<footnote n="3">³Laboratoire CEISAM-UMR CNRS 6230, Université de Nantes, 44300 Nantes, France</footnote>

<footnote>ᵃ) Author to whom correspondence should be addressed: jesper.norell@fysik.su.se</footnote>

<footnote>ᵇ) odelius@fysik.su.se</footnote>

<footnote>ᶜ) morgane.vacher@univ-nantes.fr</footnote>
# Ultrafast dynamics of photo-excited 2-thiopyridone: Theoretical insights into triplet state population and proton transfer pathways

EP

Cite as: Struct. Dyn. **7**, 024101 (2020); doi: 10.1063/1.5143228

Submitted: 20 December 2019 · Accepted: 29 February 2020 ·

Published Online: 17 March 2020

<figure><img src="image_2.png" /></figure>

View Online

<figure><img src="image_3.png" /></figure>

Export Citation

<figure><img src="image_4.png" /></figure>

CrossMark

Jesper Norell,<fnref n="1" />,ᵃ⁾ Michael Odelius,<fnref n="1" />,ᵇ⁾ and Morgane Vacher<fnref n="2" />,<fnref n="3" />,ᶜ⁾

## AFFILIATIONS

## ABSTRACT

Ultrafast non-adiabatic dynamics of the small heteroaromatic compound 2-thiopyridone has been studied with surface hopping simulations based on multi-configurational quantum chemistry. Initial excitation of the bright $S_2(\pi, \pi^*)$ state is found to promptly relax to $S_1(n, \pi^*)$ through in-plane motion. The subsequent dynamics are oppositely driven by out-of-plane motion, which results in both complex population transfers among all the available states and intersystem crossing predominantly through the “El-Sayed forbidden” $S_1(n, \pi^*)$ to $T_2(n, \pi^*)$ channel, through significant mixing of electronic excitation characters. Despite this complexity, the femto- to picosecond triplet population, expected from several spectroscopic measurements, is well described as a simple exponential decay of the singlet state manifold. No proton transfer is found in the reported trajectories, but two mechanisms for its possible mediation in previously reported experiments are proposed based on the observed structural dynamics: (i) ultrafast intra-molecular transfer driven by the initially coherent in-plane motion and (ii) inter-molecular solvent-mediated transfer driven by the out-of-plane modes that dominate the later motion.

© 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5143228

## I. INTRODUCTION

Due to the importance of light induced processes in biomolecules, in general, and in nucleic acids, in particular, the photochemistry of heteroaromatic compounds has been the subject of extensive investigations.<fnref n="1" /> Among these, thiosubstituted (organosulfuric) compounds have gathered increasing interest from ultrafast science due to their wide range of applications and rich excited-state dynamics.<fnref n="2" />⁻¹⁷ One such compound, 2-thiopyridone (2-TP, see Fig. 2), has primarily been studied as a compact model system for excited-state proton transfer (ESPT),¹⁸⁻²⁴ which is related to, for instance, photo-protection mechanisms of DNA²⁵ and pigmentation of the human skin.²⁶

While it is long established that the nitrogen protonated 2-TP can tautomerize into sulfur protonated 2-mercaptopyridine upon a change in the chemical environment (gas phase, solid, or solvated),²⁷⁻³¹ it is only recently that time-resolved measurements¹⁸,²⁰⁻²² have been

used to explicitly investigate the possibility for the corresponding UV-pumped process. In the measurement with the highest temporal resolution, traces of N-H dissociation were found below picosecond (ps) timescales using N1s resonant inelastic x-ray scattering (RIXS), but, due to limited statistics, the electronic states involved in the process could not be resolved.²⁰ Two pico- to nanosecond (ns) measurements, performed using N1s²² and S1s²¹ near-edge x-ray absorption (NEXAFS), however, instead find the ESPT to be a nanosecond process, consistent with previous vibrational measurements.¹⁸ In particular, our recent study,²² combining time-resolved N1s NEXAFS with spectrum simulations of valence-excited states, provided strong evidence for a dominant reaction pathway on these timescales, shared for $S_2$ and $S_4$ initial excitation, where ESPT results from nanosecond decay of the $T_1$ state. Notably, none of the studies have been able to address whether the ESPT occurs intra-molecularly or inter-molecularly through solvent-mediated formation of a deprotonated anion $2-TP⁻$.

<footnote n="1">¹Department of Physics, AlbaNova University Center, Stockholm University, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="2">²Department of Chemistry-Ångström Laboratory, Uppsala University, 75121 Uppsala, Sweden</footnote>

<footnote n="3">³Laboratoire CEISAM-UMR CNRS 6230, Université de Nantes, 44300 Nantes, France</footnote>

<footnote>ᵃ) Author to whom correspondence should be addressed: jesper.norell@fysik.su.se</footnote>

<footnote>ᵇ) odelius@fysik.su.se</footnote>

<footnote>ᶜ) morgane.vacher@univ-nantes.fr</footnote>
The joint body of evidence for the ESPT pathway in 2-TP currently contains several gaps, particularly with respect to the ultrafast dynamics that leads to triplet population. Highly related, as a culmination of many previous experimental<sup>2-9</sup> and theoretical<sup>10-17</sup> efforts, the groups of Crespo-Hernández and González recently offered a mechanical explanation for the general trend of rapid and efficient triplet population in thionated compounds as a stabilization that depends on the localization of the electronic excitation.<sup>42</sup> The spectroscopic measurements, based on UV absorption, were in this case supplemented with *ab initio* simulations of the non-adiabatic molecular dynamics, where recent extensions of trajectory surface hopping methods<sup>33</sup> allowed the authors to describe, with relative affordability, non-adiabatic dynamics with spin-orbit induced transitions.<sup>34-36</sup> With a computational description of the complete photo-chemical process, it can, thereby, be interpreted in well-known chemical concepts such as molecular orbitals, potential energy surfaces, and conical intersections, which may not be possible to similarly extract from measurements.

In the current work, we simulate the non-adiabatic molecular dynamics, including spin-orbit coupling induced transitions, of photo-excited 2-TP. This provides a link between the excited-state dynamics probed in previous time-resolved measurements<sup>18,20-22</sup> and quantum chemical insights gained from previous simulations.<sup>18,22-24,31</sup> Our main goal is to identify inherent aspects of the ultrafast dynamics, both electronic and structural, which result in the observed triplet population and may drive ESPT along previously reported or alluded pathways at various timescales.

## II. COMPUTATIONAL DETAILS

Non-adiabatic molecular dynamics of 2-TP was simulated by trajectory surface hopping<sup>33</sup> as implemented in the SHARC code,<sup>34,35</sup> version 1.01. The simulations included three singlet ($S_0$, $S_1$, and $S_2$) and two triplet ($T_1$ and $T_2$) states, as denoted in the basis of molecular Coulomb Hamiltonian (MCH) states (i.e., with pure spin multiplicities as obtained directly from the electronic structure calculations) with propagation performed in the spin-mixed basis of fully diagonal states. A time step of 0.5 femtoseconds (fs) was employed for the nuclear dynamics, with 100 substeps for propagation of the electronic state amplitudes. Non-adiabatic couplings were estimated from the overlap of electronic wave functions at different time steps, referred to as the local diabatic or overlap method.<sup>37-39</sup> An energy based decoherence correction<sup>40</sup> of 0.1 Hartree was further applied to correct for the artificially high coherence between different electronic states otherwise expected in trajectory surface hopping frameworks. An ensemble of 100 trajectories was simulated, with initial conditions sampled from a Wigner distribution to reproduce the vibrational ground state of $S_0$; all simulations were initialized populating the $S_2$ state.

The underlying electronic structure calculations were performed in OpenMolcas<sup>41</sup> with a CASSCF(12, 10)<sup>42</sup> setup that included: the full $\pi$ system ($\pi_1, \pi_2, \pi_3, \pi_S, \pi_1^*, \pi_2^*, \pi_3^*$) and the lone-pair $n_S$ to describe the valence excitations, together with $\sigma_{N-H}$ and $\sigma_{N-H}^*$ to describe the potential breaking of the N-H bond, as visualized in Fig. 1. The same level of theory was also applied to (i) optimize the $S_0$ geometry and calculate its vibrational modes and vertical electronic excitations and (ii) optimize pair-wise minimum energy crossing points (MECPs) [or, more specifically, minimum energy conical intersections (MECIs) for states with the same spin multiplicity] between the electronic states. Second-order perturbation theory

corrections for the vertical excitation energies in the Franck-Condon geometry were additionally calculated using multi-state CASPT2.<sup>43,44</sup> Scalar relativistic calculations using the Douglas-Kroll-Hess formulation<sup>45,46</sup> was performed employing the relativistic ANO-RCC-VDZP<sup>47</sup> basis set (which is shortly compared to the larger ANO-RCC-VTZP basis in supplementary material Sec. S1) together with Cholesky factorization.<sup>48,49</sup> Spin-orbit coupling elements were computed using an atomic mean field integral (AMFI)<sup>50</sup> approach within the RASSI framework.<sup>51,52</sup> Molecular orbitals were visualized with the luscus program.<sup>53</sup>

Further details of the complete computational setup, in the form of example inputs, are provided in the supplementary material.

## III. RESULTS AND DISCUSSION

### A. Electronic states

To describe and rationalize the non-adiabatic dynamics, we first introduce the electronic states that will be involved. An overview of the first four excited states of 2-TP and their potential energy surfaces in the region of the explored dynamics, schematically constructed in accordance with our later drawn conclusions, are given in Fig. 2. The electronic excitation energies computed in the Franck-Condon geometry at CASSCF and CASPT2 levels of theory are gathered in Table I(a). The singlet manifold consists of the dark $S_1(n, \pi^*)$ and the bright $S_2(\pi, \pi^*)$, whereas the triplet manifold contains $T_1(\pi, \pi^*)$ and $T_2(n, \pi^*)$. Including dynamic correlation with perturbation theory leads to energy differences of less than 0.2 eV for singlet states and a bit larger for triplet states. Importantly, the electronic state characters and their energy ordering remain the same upon adding dynamic correlation. It is interesting to remark that the energy ordering of excitation characters $((n, \pi^*)$ and $(\pi, \pi^*))$ is inverted between the two manifolds (singlet and triplet) and all states are found in a narrow energy range, in particular, the lowest three, at both levels of theory. While the excitation characters may be strictly identified for all states in the currently presented Franck-Condon analysis, we note that the same does not hold through the later presented trajectories. For discussion of the non-adiabatic dynamics, the states will, therefore, in general, simply be denoted by their energy ordered MCH state labels, i.e., without specified excitation character.

In a single-electron picture, all the states nominally result from excitation into the lowest unoccupied molecular orbital (LUMO) $\pi_1^*$,

<figure><img src="image_4.png" /><figcaption>FIG. 1. Active molecular orbitals in the CASSCF calculations, here shown for the $S_0$ optimized geometry. The studied valence states can nominally be assigned to excitation from the sulfur dominated $\pi_S$ ($S_2$ and $T_1$) and the sulfur lone-pair $n_S$ ($S_1$ and $T_2$) into $\pi_1^*$.</figcaption></figure>
The joint body of evidence for the ESPT pathway in 2-TP currently contains several gaps, particularly with respect to the ultrafast dynamics that leads to triplet population. Highly related, as a culmination of many previous experimental<sup>2-9</sup> and theoretical<sup>10-17</sup> efforts, the groups of Crespo-Hernández and González recently offered a mechanistic explanation for the general trend of rapid and efficient triplet population in thionated compounds as a stabilization that depends on the localization of the electronic excitation.<sup>32</sup> The spectroscopic measurements, based on UV absorption, were in this case supplemented with *ab initio* simulations of the non-adiabatic molecular dynamics, where recent extensions of trajectory surface hopping methods<sup>33</sup> allowed the authors to describe, with relative affordability, non-adiabatic dynamics with spin-orbit induced transitions.<sup>34-36</sup> With a computational description of the complete photo-chemical process, it can, thereby, be interpreted in well-known chemical concepts such as molecular orbitals, potential energy surfaces, and conical intersections, which may not be possible to similarly extract from measurements.

In the current work, we simulate the non-adiabatic molecular dynamics, including spin-orbit coupling induced transitions, of photo-excited 2-TP. This provides a link between the excited-state dynamics probed in previous time-resolved measurements<sup>18,20-22</sup> and quantum chemical insights gained from previous simulations.<sup>18,22-24,31</sup> Our main goal is to identify inherent aspects of the ultrafast dynamics, both electronic and structural, which result in the observed triplet population and may drive ESPT along previously reported or alluded pathways at various timescales.

## II. COMPUTATIONAL DETAILS

Non-adiabatic molecular dynamics of 2-TP was simulated by trajectory surface hopping<sup>33</sup> as implemented in the SHARC code,<sup>34,35</sup> version 1.01. The simulations included three singlet ($S_0$, $S_1$, and $S_2$) and two triplet ($T_1$ and $T_2$) states, as denoted in the basis of molecular Coulomb Hamiltonian (MCH) states (i.e., with pure spin multiplicities as obtained directly from the electronic structure calculations) with propagation performed in the spin-mixed basis of fully diagonal states. A time step of 0.5 femtoseconds (fs) was employed for the nuclear dynamics, with 100 substeps for propagation of the electronic state amplitudes. Non-adiabatic couplings were estimated from the overlap of electronic wave functions at different time steps, referred to as the local diabatic or overlap method.<sup>37-39</sup> An energy based decoherence correction<sup>40</sup> of 0.1 Hartree was further applied to correct for the artificially high coherence between different electronic states otherwise expected in trajectory surface hopping frameworks. An ensemble of 100 trajectories was simulated, with initial conditions sampled from a Wigner distribution to reproduce the vibrational ground state of $S_0$; all simulations were initialized populating the $S_2$ state.

The underlying electronic structure calculations were performed in OpenMolcas<sup>41</sup> with a CASSCF(12, 10)<sup>42</sup> setup that included: the full $\pi$ system ($\pi_1, \pi_2, \pi_3, \pi_S, \pi_1^*, \pi_2^*, \pi_3^*$) and the lone-pair $n_S$ to describe the valence excitations, together with $\sigma_{N-H}$ and $\sigma_{N-H}^*$ to describe the potential breaking of the N-H bond, as visualized in Fig. 1. The same level of theory was also applied to (i) optimize the $S_0$ geometry and calculate its vibrational modes and vertical electronic excitations and (ii) optimize pair-wise minimum energy crossing points (MECPs) [or, more specifically, minimum energy conical intersections (MECIs) for states with the same spin multiplicity] between the electronic states. Second-order perturbation theory

corrections for the vertical excitation energies in the Franck-Condon geometry were additionally calculated using multi-state CASPT2.<sup>43,44</sup> Scalar relativistic calculations using the Douglas-Kroll-Hess formulation<sup>45,46</sup> was performed employing the relativistic ANO-RCC-VDZP<sup>47</sup> basis set (which is shortly compared to the larger ANO-RCC-VTZP basis in supplementary material Sec. S1) together with Cholesky factorization.<sup>48,49</sup> Spin-orbit coupling elements were computed using an atomic mean field integral (AMFI)<sup>50</sup> approach within the RASSI framework.<sup>51,52</sup> Molecular orbitals were visualized with the luscus program.<sup>53</sup>

Further details of the complete computational setup, in the form of example inputs, are provided in the supplementary material.

## III. RESULTS AND DISCUSSION

### A. Electronic states

To describe and rationalize the non-adiabatic dynamics, we first introduce the electronic states that will be involved. An overview of the first four excited states of 2-TP and their potential energy surfaces in the region of the explored dynamics, schematically constructed in accordance with our later drawn conclusions, are given in Fig. 2. The electronic excitation energies computed in the Franck-Condon geometry at CASSCF and CASPT2 levels of theory are gathered in Table I(a). The singlet manifold consists of the dark $S_1(n, \pi^*)$ and the bright $S_2(\pi, \pi^*)$, whereas the triplet manifold contains $T_1(\pi, \pi^*)$ and $T_2(n, \pi^*)$. Including dynamic correlation with perturbation theory leads to energy differences of less than 0.2 eV for singlet states and a bit larger for triplet states. Importantly, the electronic state characters and their energy ordering remain the same upon adding dynamic correlation. It is interesting to remark that the energy ordering of excitation characters ((n, π^*) and (π, π^*)) is inverted between the two manifolds (singlet and triplet) and all states are found in a narrow energy range, in particular, the lowest three, at both levels of theory. While the excitation characters may be strictly identified for all states in the currently presented Franck-Condon analysis, we note that the same does not hold through the later presented trajectories. For discussion of the non-adiabatic dynamics, the states will, therefore, in general, simply be denoted by their energy ordered MCH state labels, i.e., without specified excitation character.

In a single-electron picture, all the states nominally result from excitation into the lowest unoccupied molecular orbital (LUMO) $\pi_1^*$,

<figure><img src="image_5.png" /><figcaption>FIG. 1. Active molecular orbitals in the CASSCF calculations, here shown for the $S_0$ optimized geometry. The studied valence states can nominally be assigned to excitation from the sulfur dominated $\pi_S$ ($S_2$ and $T_1$) and the sulfur lone-pair $n_S$ ($S_1$ and $T_2$) into $\pi_1^*$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Kinetic model schematic</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIG. 2.** Schematic potential energy surfaces in the region of 2-TP photo-excitation (left). Schematic of the kinetic model derived from the simulations and fitted to the electronic state populations. Solid lines indicate the main pathway and dashed lines the minority pathways (right).</figcaption></figure>

$S_2$ and $T_1$ from the sulfur-dominated $\pi_S$, $S_1$ and $T_2$ from the in-plane sulfur lone-pair $n_S$. Of particular relevance for later interpretation, we note that $\pi_S$ and $n_S$ correspond to the nearly degenerate highest occupied molecular orbital (HOMO) and HOMO-1 in a single-determinant description, respectively,²³ where the sulfur dominance of the former makes it comparable to an out-of-plane lone-pair. Visualizations of all the active molecular orbitals are found in Fig. 1.

## B. Electronic dynamics

Non-adiabatic molecular dynamics simulations were performed to investigate the relaxation upon the population of the $S_2$ bright state, with CASSCF as the underlying electronic structure method. A number of studies, based on *ab initio* multiple spawning,⁵⁴⁻⁶³ in particular, but also on other methods,¹⁷,⁶⁴ have previously demonstrated the ability of CASPT2 based simulations of non-adiabatic dynamics to yield

<table><caption>**TABLE 1.** (a) Excitation energies of the electronic states in the Franck-Condon region. Upper row: directly from CASSCF calculations, as employed for the dynamics. Lower row in parentheses: with added multistate CASPT2 correction. (b) Total number of surface hops from the old MCH state (rows) to a new MCH state (columns) found in the 100 trajectories. Only hops with 10 or more occurrences are shown, and the forward transitions of the main pathway are emphasized in bold.</caption><caption>(a) Excitation energies in units of eV.</caption><thead><tr><th>S<sub>1</sub>(n, π<sup>*</sup>)</th><th>S<sub>2</sub>(π, π<sup>*</sup>)</th><th>T<sub>1</sub>(π, π<sup>*</sup>)</th><th>T<sub>2</sub>(n, π<sup>*</sup>)</th></tr></thead><tbody><tr><td>3.16<br>(3.34)</td><td>3.82<br>(3.64)</td><td>2.64<br>(2.92)</td><td>2.90<br>(3.31)</td></tr></tbody></table>

<table><caption>(b) Total number of MCH surface hops.</caption><thead><tr><th>From To</th><th>S<sub>0</sub></th><th>S<sub>1</sub></th><th>S<sub>2</sub></th><th>T<sub>1</sub></th><th>T<sub>2</sub></th></tr></thead><tbody><tr><th>S<sub>0</sub></th><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><th>S<sub>1</sub></th><td>13</td><td>...</td><td>348</td><td>...</td><td>298</td></tr><tr><th>S<sub>2</sub></th><td>...</td><td>441</td><td>...</td><td>...</td><td>10</td></tr><tr><th>T<sub>1</sub></th><td>...</td><td>...</td><td>...</td><td>...</td><td>4849</td></tr><tr><th>T<sub>2</sub></th><td>...</td><td>251</td><td>...</td><td>4869</td><td>...</td></tr></tbody></table>

even quantitative results. Yet, due to both the difficulty in the implementation of CASPT2 gradients and their steep computational cost, simulations that include spin-orbit coupling have instead mainly been based on the CASSCF level of theory.⁶⁵⁻⁶⁸ Due to the lack of dynamic correlation in the current description, later conclusions will primarily be drawn about qualitative rather than quantitative aspects of the simulated dynamics. The simulated trajectories were all started in the bright $S_2(\pi, \pi^*)$ state, which has previously been experimentally pumped with 343 nm (Ref. 22) and 400 nm (Ref. 21) lasers and shares a common reaction pathway on picosecond timescales and longer with the second bright $S_4(\pi, \pi^*)$ state,²² pumped instead with 258 nm (Ref. 18) and 266 nm (Ref. 22) lasers.

From the ensemble of trajectories, we first consider the kinetics of the electronic state population transfers. The relative populations of the MCH states as a function of time are shown in Fig. 3 together with a kinetic model fit. Comparison of different state representations (see supplementary material Sec. S3) displays only small differences for the overall evolution of the populations, which thereby demonstrates a stability of the populations themselves despite later discussed complications in the electronic state characters. The main clues to the reaction pathway are found in Table I(b), which enumerates the number of surface hops between the states. The dominant pathway is, in this way, readily identified as a sequential decay in reverse energy order, i.e., $S_2$ to $S_1$ to $T_2$ to $T_1$, hereafter referred to as the “main” pathway. Yet, the total kinetics cannot be so easily described for at least two reasons: first, because of the numerous “backward” transitions within the main pathway, which disallows the population transfer to be accurately modeled with a simple exponential “forward” decay and second, due to the existence of several “minority” pathways, which, for instance, partially re-populates $S_0$ at very early times.

Our attempts to fit the obtained populations with a kinetic model are described in more detail in the supplementary material, Sec. S4. In short, we find that a model with nine different rate constants (three main forward, three main backward, and three minority forward transitions), as schematically depicted in Fig. 2 and with the resulting fit shown in Fig. 3, is necessary to capture all qualitative aspects of the kinetics. Consequently, the lifetime of individual electronic states can

<figure><img src="image_7.png" /><figcaption>**FIG. 3.** Electronic state populations from the 100 trajectory ensemble (thick shaded lines) compared to the fit obtained with the full kinetic model (thin opaque lines, see supplementary material Sec. S4 for further details).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Kinetic model schematic</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**FIG. 2.** Schematic potential energy surfaces in the region of 2-TP photo-excitation (left). Schematic of the kinetic model derived from the simulations and fitted to the electronic state populations. Solid lines indicate the main pathway and dashed lines the minority pathways (right).</figcaption></figure>

$S_2$ and $T_1$ from the sulfur-dominated $\pi_S$, $S_1$ and $T_2$ from the in-plane sulfur lone-pair $n_S$. Of particular relevance for later interpretation, we note that $\pi_S$ and $n_S$ correspond to the nearly degenerate highest occupied molecular orbital (HOMO) and HOMO-1 in a single-determinant description, respectively,²³ where the sulfur dominance of the former makes it comparable to an out-of-plane lone-pair. Visualizations of all the active molecular orbitals are found in Fig. 1.

## B. Electronic dynamics

Non-adiabatic molecular dynamics simulations were performed to investigate the relaxation upon the population of the $S_2$ bright state, with CASSCF as the underlying electronic structure method. A number of studies, based on *ab initio* multiple spawning,⁵⁴⁻⁶³ in particular, but also on other methods,¹⁷,⁶⁴ have previously demonstrated the ability of CASPT2 based simulations of non-adiabatic dynamics to yield

<table><caption>**TABLE 1.** (a) Excitation energies of the electronic states in the Franck-Condon region. Upper row: directly from CASSCF calculations, as employed for the dynamics. Lower row in parentheses: with added multistate CASPT2 correction. (b) Total number of surface hops from the old MCH state (rows) to a new MCH state (columns) found in the 100 trajectories. Only hops with 10 or more occurrences are shown, and the forward transitions of the main pathway are emphasized in bold.</caption><caption>(a) Excitation energies in units of eV.</caption><thead><tr><th>S<sub>1</sub>(n, π<sup>*</sup>)</th><th>S<sub>2</sub>(π, π<sup>*</sup>)</th><th>T<sub>1</sub>(π, π<sup>*</sup>)</th><th>T<sub>2</sub>(n, π<sup>*</sup>)</th></tr></thead><tbody><tr><td>3.16<br>(3.34)</td><td>3.82<br>(3.64)</td><td>2.64<br>(2.92)</td><td>2.90<br>(3.31)</td></tr></tbody></table>

<table><caption>(b) Total number of MCH surface hops.</caption><thead><tr><th>From To</th><th>S<sub>0</sub></th><th>S<sub>1</sub></th><th>S<sub>2</sub></th><th>T<sub>1</sub></th><th>T<sub>2</sub></th></tr></thead><tbody><tr><td>S<sub>0</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>S<sub>1</sub></td><td>13</td><td>...</td><td>348</td><td>...</td><td>298</td></tr><tr><td>S<sub>2</sub></td><td>...</td><td>441</td><td>...</td><td>...</td><td>10</td></tr><tr><td>T<sub>1</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>4849</td></tr><tr><td>T<sub>2</sub></td><td>...</td><td>251</td><td>...</td><td>4869</td><td>...</td></tr></tbody></table>

even quantitative results. Yet, due to both the difficulty in the implementation of CASPT2 gradients and their steep computational cost, simulations that include spin-orbit coupling have instead mainly been based on the CASSCF level of theory.⁶⁵⁻⁶⁸ Due to the lack of dynamic correlation in the current description, later conclusions will primarily be drawn about qualitative rather than quantitative aspects of the simulated dynamics. The simulated trajectories were all started in the bright $S_2(\pi, \pi^*)$ state, which has previously been experimentally pumped with 343 nm (Ref. 22) and 400 nm (Ref. 21) lasers and shares a common reaction pathway on picosecond timescales and longer with the second bright $S_4(\pi, \pi^*)$ state,²² pumped instead with 258 nm (Ref. 18) and 266 nm (Ref. 22) lasers.

From the ensemble of trajectories, we first consider the kinetics of the electronic state population transfers. The relative populations of the MCH states as a function of time are shown in Fig. 3 together with a kinetic model fit. Comparison of different state representations (see supplementary material Sec. S3) displays only small differences for the overall evolution of the populations, which thereby demonstrates a stability of the populations themselves despite later discussed complications in the electronic state characters. The main clues to the reaction pathway are found in Table I(b), which enumerates the number of surface hops between the states. The dominant pathway is, in this way, readily identified as a sequential decay in reverse energy order, i.e., $S_2$ to $S_1$ to $T_2$ to $T_1$, hereafter referred to as the “main” pathway. Yet, the total kinetics cannot be so easily described for at least two reasons: first, because of the numerous “backward” transitions within the main pathway, which disallows the population transfer to be accurately modeled with a simple exponential “forward” decay and second, due to the existence of several “minority” pathways, which, for instance, partially re-populates $S_0$ at very early times.

Our attempts to fit the obtained populations with a kinetic model are described in more detail in the supplementary material, Sec. S4. In short, we find that a model with nine different rate constants (three main forward, three main backward, and three minority forward transitions), as schematically depicted in Fig. 2 and with the resulting fit shown in Fig. 3, is necessary to capture all qualitative aspects of the kinetics. Consequently, the lifetime of individual electronic states can

<figure><img src="image_8.png" /><figcaption>**FIG. 3.** Electronic state populations from the 100 trajectory ensemble (thick shaded lines) compared to the fit obtained with the full kinetic model (thin opaque lines, see supplementary material Sec. S4 for further details).</figcaption></figure>
be neither simply nor robustly extracted, as they result from a complex balance between backward and forward transfers of population and are highly sensitive to changes of the fitted model. We assign the complexity in the population kinetics to primarily result from large similarities between the states, both in their excitation energies and excitation characters, as elaborated in this section and Sec. III C.

The small energy splitting of the electronic states, already noted from their vertical excitations, is presumably the main cause for the significant contribution of backward transitions within the main pathway. An estimate of how the energy splittings develop over time, based on the spin-orbit coupled state representation wherein the SHARC dynamics are propagated, is presented in supplementary material Sec. S2. In short, following large fluctuations within the first ~100 fs, the mean values of the splittings stabilize toward roughly 0.4 eV ($S_1$ to $T_2$) and 0.3 eV ($T_1$ to $T_2$). As a result, the $S_1$, $T_1$, and $T_2$ states share comparable populations at 1 ps, in contrast to multiple studies of similar systems where a single $T_1$ state dominates at such timescales.¹⁴,¹⁷,³² The increased tendency for backward transitions can also arise from limited “funneling” toward the associated conical intersection in the “upper” state and subsequently away from it in the “lower” state, as previously observed for $LiH_2$.⁶⁹ Indeed, comparing the MECIs (which represents well the hopping geometries as shown in Sec. III C) for $S_2-S_1$ and $T_2-T_1$, the former is categorized as peaked and bifurcating and the latter as a sloped single-path.⁷⁰ That is, the singlet transition (with fewer backward occurrences) takes place around a minimum of the upper state surface with multiple preferred relaxation directions on the lower surface (which should both promote funneling), whereas the triplet transition (with more backward occurrences) takes place on a slope of the upper state with a single downwards pathway in the lower state (which may both inhibit funneling). The tendency for backward transfer of population, in general, would likely be reduced in a solvent by rapid dissipation of kinetic energy and thereby not be directly relatable to previous experiments.¹⁸,²⁰⁻²²

In contrast to the energies, the similarities in excitation character cannot at all be captured within a Franck-Condon picture. The effect is instead most clearly demonstrated by how intersystem crossing (ISC) is in the current case dominated by $S_1$ to $T_2$ transfer, which is El-Sayed forbidden as the two states share the same nominal electronic configuration.⁷¹,⁷² The ISC must, therefore, be driven by symmetry-breaking in the structural dynamics, where the already noted similarity in both energy and localization of the $n_S$ and $\pi_S$ orbitals then allows for significant mixing of the electronic characters. Based on a single point calculation of the MECP geometry between the $S_1$ and $T_2$ states (which represents well the hopping geometries as shown in Sec. III C), we find as leading terms of the configuration expansions for the states $S_1 = 0.69(n_S, \pi^+) + 0.58(\pi_S, \pi^+)$ and $T_2 = 0.73(n_S, \pi^+) - 0.58(\pi_S, \pi^+)$. The two states are, thereby, essentially orthogonal in their spatial wave functions, similarly as expected for the “El-Sayed allowed” transitions, but consist, in contrast, of heavily mixed contributions from in-plane and out-of-plane excitation characters. Similarly, for a sampling over five randomly selected hop-inducing geometries (as introduced in the Sec. III C) for the same transition, the weight of the dominant configuration has a mean (standard deviation) of 0.63 (±0.13) and 0.81 (±0.14) for the $S_1$ and $T_2$ states, respectively. It can, for these reasons, also be assumed that all the electronic states, as they all nominally involve excitations between the same orbitals, similarly carry a mixed character throughout parts of the non-adiabatic dynamics.

While the kinetics of the individual electronic states appear to be difficult to model, a simple trend emerges when they are summed together according to their spin multiplicities. The decay of the singlet manifold population, as shown in Fig. 4, can be almost as well described with a single exponential fit as with the full kinetic model. The simplified fit allows us to better quantify the triplet state population transfer, with the rate constant that translates into an ~1 ps lifetime of the singlet manifold. The result is consistent with x-ray measurements in solutions, where $T_1$ appears to be dominantly populated within the time resolution of ~20 ps.²¹,²² The trend of the efficient triplet population upon thionation observed for a number of compounds,³² thereby, appears to hold also with the active involvement of more electronic states in the reaction pathway. Further extrapolated within the 2-TP system, this trend would explain why the same reaction pathway was also observed for $S_4$ excitation at pico- to nanosecond scales,²² through similar ultrafast relaxation to the triplet manifold.

### C. Structural dynamics

With an overview of the population transfer kinetics, we now turn to the structural dynamics. To rationalize how the structural distortions drive the electronic transitions, we first consider the “hop-inducing” geometries. Figure 5 shows, for each forward transition within the main and the minority pathways, the mean value $\mu$ and standard deviation $\sigma$ of all geometries that precede the corresponding surface hops. The backward transitions, as shown in supplementary material Fig. S6, show no appreciable differences. The geometries are expressed in normal mode coordinates, some of which are visualized in Fig. 6, while the rest can be found in supplementary material Figs. S7–S9. The data are further compared with pair-wise MECP geometries, visualized in supplementary material Fig. S5.

Starting with the main pathway, the initial internal conversion (IC) from $S_2$ to $S_1$ occurs through minimal structural distortion. This can be rationalized by their small energy gap already in the Franck-Condon region and is consistent with the ultrafast population transfer from $S_2$ to $S_1$ (see Fig. 3). In contrast, the $S_1$ to $T_2$ ISC is driven

<figure><img src="image_8.png" /><figcaption>**FIG. 4.** Singlet and triplet state populations (summed over individual MCH states, thick shaded lines) compared to the full kinetic model fit (thin full lines, see supplementary material Sec. S4 for further details) and a simple exponential fit with a lifetime of $\tau = 1058$ fs ~1 ps of the singlet manifold (dashed lines).</figcaption></figure>
be neither simply nor robustly extracted, as they result from a complex balance between backward and forward transfers of population and are highly sensitive to changes of the fitted model. We assign the complexity in the population kinetics to primarily result from large similarities between the states, both in their excitation energies and excitation characters, as elaborated in this section and Sec. III C.

The small energy splitting of the electronic states, already noted from their vertical excitations, is presumably the main cause for the significant contribution of backward transitions within the main pathway. An estimate of how the energy splittings develop over time, based on the spin-orbit coupled state representation wherein the SHARC dynamics are propagated, is presented in supplementary material Sec. S2. In short, following large fluctuations within the first ~100 fs, the mean values of the splittings stabilize toward roughly 0.4 eV ($S_1$ to $T_2$) and 0.3 eV ($T_1$ to $T_2$). As a result, the $S_1$, $T_1$, and $T_2$ states share comparable populations at 1 ps, in contrast to multiple studies of similar systems where a single $T_1$ state dominates at such timescales.¹⁴,¹⁷,³² The increased tendency for backward transitions can also arise from limited “funneling” toward the associated conical intersection in the “upper” state and subsequently away from it in the “lower” state, as previously observed for $LiH_2$.⁶⁹ Indeed, comparing the MECIs (which represents well the hopping geometries as shown in Sec. III C) for $S_2-S_1$ and $T_2-T_1$, the former is categorized as peaked and bifurcating and the latter as a sloped single-path.⁷⁰ That is, the singlet transition (with fewer backward occurrences) takes place around a minimum of the upper state surface with multiple preferred relaxation directions on the lower surface (which should both promote funneling), whereas the triplet transition (with more backward occurrences) takes place on a slope of the upper state with a single downwards pathway in the lower state (which may both inhibit funneling). The tendency for backward transfer of population, in general, would likely be reduced in a solvent by rapid dissipation of kinetic energy and thereby not be directly relatable to previous experiments.¹⁸,²⁰⁻²²

In contrast to the energies, the similarities in excitation character cannot at all be captured within a Franck-Condon picture. The effect is instead most clearly demonstrated by how intersystem crossing (ISC) is in the current case dominated by $S_1$ to $T_2$ transfer, which is El-Sayed forbidden as the two states share the same nominal electronic configuration.⁷¹,⁷² The ISC must, therefore, be driven by symmetry-breaking in the structural dynamics, where the already noted similarity in both energy and localization of the $n_S$ and $\pi_S$ orbitals then allows for significant mixing of the electronic characters. Based on a single point calculation of the MECP geometry between the $S_1$ and $T_2$ states (which represents well the hopping geometries as shown in Sec. III C), we find as leading terms of the configuration expansions for the states $S_1 = 0.69(n_S, \pi^+) + 0.58(\pi_S, \pi^+)$ and $T_2 = 0.73(n_S, \pi^+) - 0.58(\pi_S, \pi^+)$. The two states are, thereby, essentially orthogonal in their spatial wave functions, similarly as expected for the “El-Sayed allowed” transitions, but consist, in contrast, of heavily mixed contributions from in-plane and out-of-plane excitation characters. Similarly, for a sampling over five randomly selected hop-inducing geometries (as introduced in the Sec. III C) for the same transition, the weight of the dominant configuration has a mean (standard deviation) of 0.63 (±0.13) and 0.81 (±0.14) for the $S_1$ and $T_2$ states, respectively. It can, for these reasons, also be assumed that all the electronic states, as they all nominally involve excitations between the same orbitals, similarly carry a mixed character throughout parts of the non-adiabatic dynamics.

While the kinetics of the individual electronic states appear to be difficult to model, a simple trend emerges when they are summed together according to their spin multiplicities. The decay of the singlet manifold population, as shown in Fig. 4, can be almost as well described with a single exponential fit as with the full kinetic model. The simplified fit allows us to better quantify the triplet state population transfer, with the rate constant that translates into an ~1 ps lifetime of the singlet manifold. The result is consistent with x-ray measurements in solutions, where $T_1$ appears to be dominantly populated within the time resolution of ~20 ps.²¹,²² The trend of the efficient triplet population upon thionation observed for a number of compounds,³² thereby, appears to hold also with the active involvement of more electronic states in the reaction pathway. Further extrapolated within the 2-TP system, this trend would explain why the same reaction pathway was also observed for $S_4$ excitation at pico- to nanosecond scales,²² through similar ultrafast relaxation to the triplet manifold.

### C. Structural dynamics

With an overview of the population transfer kinetics, we now turn to the structural dynamics. To rationalize how the structural distortions drive the electronic transitions, we first consider the “hop-inducing” geometries. Figure 5 shows, for each forward transition within the main and the minority pathways, the mean value $\mu$ and standard deviation $\sigma$ of all geometries that precede the corresponding surface hops. The backward transitions, as shown in supplementary material Fig. S6, show no appreciable differences. The geometries are expressed in normal mode coordinates, some of which are visualized in Fig. 6, while the rest can be found in supplementary material Figs. S7–S9. The data are further compared with pair-wise MECP geometries, visualized in supplementary material Fig. S5.

Starting with the main pathway, the initial internal conversion (IC) from $S_2$ to $S_1$ occurs through minimal structural distortion. This can be rationalized by their small energy gap already in the Franck-Condon region and is consistent with the ultrafast population transfer from $S_2$ to $S_1$ (see Fig. 3). In contrast, the $S_1$ to $T_2$ ISC is driven

<figure><img src="image_9.png" /><figcaption>**FIG. 4.** Singlet and triplet state populations (summed over individual MCH states, thick shaded lines) compared to the full kinetic model fit (thin full lines, see supplementary material Sec. S4 for further details) and a simple exponential fit with a lifetime of $\tau = 1058$ fs ~1 ps of the singlet manifold (dashed lines).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 5. Transition-inducing molecular geometries: comparison of the molecular geometries that precede a surface hop in the SHARC simulations (hop) to MECP geometries. The geometries are expressed as distortions along the normal mode coordinates of the electronic ground state, with each data point as the projection (12 × 3 dimensional scalar product) $p_i = \langle \Delta\vec{r}_{geo} \cdot \Delta\vec{r}_{NM,i} \rangle$ between the molecular geometry distortion $\Delta\vec{r}_{geo} = \vec{r}_{geo} - \vec{r}_{Franck-Condons}$ and the normalized normal mode displacement vector $\Delta\vec{r}_{NM,i}$ of mode i. Green dots mark in-plane modes and blue dots out-of-plane modes. The out-of-plane modes were, for each geometry, if needed, aligned to a positive phase of mode 8 by a mirroring through the molecular plane, to facilitate comparison to the MECPs.</figcaption></figure>

by strong distortions in almost all out-of-plane modes (blue dots in Fig. 5). This is, as previously discussed, necessary to enable the spin-flip through mixing of electronic state characters, which would, otherwise, be highly inefficient between a singlet and triplet state of the same configuration.⁷¹,⁷² A large number of apparent transitions within the triplet manifold then proceed through similar out-of-plane motion. Yet, considering the minuscule energy separation of the T₁ and T₂ states already in the Franck–Condon region and the character mixing induced by the symmetry breaking, it is likely that most of these reflect a rapid and erratic energy re-ordering rather than transitions between well distinguishable states.

The first minority pathway results from S₂ to T₂ ISC through small structural distortions. It, thereby, represents the “El-Sayed efficient”⁷¹,⁷² (π, π*) to (n, π*) ISC, in contrast to the El-Sayed forbidden (n, π*) to (n, π*) ISC of the main pathway. In further contrast, this ISC channel also competes with the S₂ to S₁ IC, as S₂ is not the lowest state within the manifold of excited singlets. An appreciable ISC yield already from the S₂ state, thereby, appears to directly reflect the increased ISC rate by the heavy atom effect of the sulfur. Still, a significant majority of ISC occur after relaxation within the singlet manifold to S₁, as traditionally expected.

The second minority pathway results from S₁ to S₀ IC. To separate this from the main pathway that instead proceeds through ISC, we denote the minority pathway trajectories with an alternative species S'₁. While the split between S₁ and S'₁ is for fitting purposes implemented in the decay from S₂, we deduce that it actually occurs in the subsequent dynamics on the S₁ surface as presented in Fig. 2. This as S₁ and S'₁ are reached through highly comparable geometries [compare Figs. 5(a) and 5(d)] but decays through more distinguishable geometries [compare Figs. 5(b) and 5(e)] on different timescales: S'₁ to S₀ occurs on average after 114 fs, whereas S₁ to T₂ occurs on average after

369 fs. The splitting is, therefore, likely to be (partially) induced by the bifurcating character of the S₂–S₁ MECI, mentioned in Sec. III B.

The dynamical picture based on the hop-inducing geometries can be compared with a static perspective provided by the MECP geometries also presented in Fig. 5. We note that the MECP data for out-of-plane modes, due to the planar symmetry that makes their sign arbitrary, must be compared to the total distortions enabled by both the mean values and standard deviations of the hop statistics. The exception to this is mode 8 and any other mode coherently phased to it, to which the signs were aligned to facilitate comparison. For the main pathway, we mainly note that the results appear to be in overall acceptable agreement, i.e., the transitions appear to, in most cases, occur “close” to the MECP, which thereby validates the argumentation about electronic state characters for the associated ICs and ISC invoked in the last section, based on single point calculations at the MECPs (MECIs). Regrettably, the S₂ to T₂ MECP of the first minority pathway could not be converged, due to the complications that arise from S₁ being energetically situated between the two states. Finally, the S₁ to S₀ MECI predicts much larger distortion, in several modes, than seen in the dynamics. This can be understood by how the transitions, on average, occur for a potential energy of 1.33 eV above the MECI and results in a 0.69 eV increase in kinetic energy, which shows that the transitions typically occur far from both the MECI and any degeneracy between the two surfaces.

Finally, we consider the structural dynamics also in a fully time-dependent picture. Figure 6 presents the time-evolution of selected normal mode coordinates (i.e., those that display notable changes in amplitude over time and/or clear coherence) over the whole ensemble of trajectories. To gain some measure of the coherence between trajectories, we compare the mean values of the signed and unsigned coordinates. Equal amplitudes (of, e.g., a periodic oscillation or systematic
<figure><img src="image_10.png" /><figcaption>**FIG. 5.** Transition-inducing molecular geometries: comparison of the molecular geometries that precede a surface hop in the SHARC simulations (hop) to MECP geometries. The geometries are expressed as distortions along the normal mode coordinates of the electronic ground state, with each data point as the projection (12 × 3 dimensional scalar product) $p_i = \langle \Delta\vec{r}_{geo} \cdot \Delta\vec{r}_{NM,i} \rangle$ between the molecular geometry distortion $\Delta\vec{r}_{geo} = \vec{r}_{geo} - \vec{r}_{Franck-Condons}$ and the normalized normal mode displacement vector $\Delta\vec{r}_{NM,i}$ of mode *i*. Green dots mark in-plane modes and blue dots out-of-plane modes. The out-of-plane modes were, for each geometry, if needed, aligned to a positive phase of mode 8 by a mirroring through the molecular plane, to facilitate comparison to the MECPs.</figcaption></figure>

by strong distortions in almost all out-of-plane modes (blue dots in Fig. 5). This is, as previously discussed, necessary to enable the spin-flip through mixing of electronic state characters, which would, otherwise, be highly inefficient between a singlet and triplet state of the same configuration.⁷¹,⁷² A large number of apparent transitions within the triplet manifold then proceed through similar out-of-plane motion. Yet, considering the minuscule energy separation of the T₁ and T₂ states already in the Franck–Condon region and the character mixing induced by the symmetry breaking, it is likely that most of these reflect a rapid and erratic energy re-ordering rather than transitions between well distinguishable states.

The first minority pathway results from S₂ to T₂ ISC through small structural distortions. It, thereby, represents the “El-Sayed efficient”⁷¹,⁷² (π, π*) to (n, π*) ISC, in contrast to the El-Sayed forbidden (n, π*) to (n, π*) ISC of the main pathway. In further contrast, this ISC channel also competes with the S₂ to S₁ IC, as S₂ is not the lowest state within the manifold of excited singlets. An appreciable ISC yield already from the S₂ state, thereby, appears to directly reflect the increased ISC rate by the heavy atom effect of the sulfur. Still, a significant majority of ISC occur after relaxation within the singlet manifold to S₁, as traditionally expected.

The second minority pathway results from S₁ to S₀ IC. To separate this from the main pathway that instead proceeds through ISC, we denote the minority pathway trajectories with an alternative species S'₁. While the split between S₁ and S'₁ is for fitting purposes implemented in the decay from S₂, we deduce that it actually occurs in the subsequent dynamics on the S₁ surface as presented in Fig. 2. This as S₁ and S'₁ are reached through highly comparable geometries [compare Figs. 5(a) and 5(d)] but decays through more distinguishable geometries [compare Figs. 5(b) and 5(e)] on different timescales: S'₁ to S₀ occurs on average after 114 fs, whereas S₁ to T₂ occurs on average after

369 fs. The splitting is, therefore, likely to be (partially) induced by the bifurcating character of the S₂–S₁ MECI, mentioned in Sec. III B.

The dynamical picture based on the hop-inducing geometries can be compared with a static perspective provided by the MECP geometries also presented in Fig. 5. We note that the MECP data for out-of-plane modes, due to the planar symmetry that makes their sign arbitrary, must be compared to the total distortions enabled by both the mean values and standard deviations of the hop statistics. The exception to this is mode 8 and any other mode coherently phased to it, to which the signs were aligned to facilitate comparison. For the main pathway, we mainly note that the results appear to be in overall acceptable agreement, i.e., the transitions appear to, in most cases, occur “close” to the MECP, which thereby validates the argumentation about electronic state characters for the associated ICs and ISC invoked in the last section, based on single point calculations at the MECPs (MECIs). Regrettably, the S₂ to T₂ MECP of the first minority pathway could not be converged, due to the complications that arise from S₁ being energetically situated between the two states. Finally, the S₁ to S₀ MECI predicts much larger distortion, in several modes, than seen in the dynamics. This can be understood by how the transitions, on average, occur for a potential energy of 1.33 eV above the MECI and results in a 0.69 eV increase in kinetic energy, which shows that the transitions typically occur far from both the MECI and any degeneracy between the two surfaces.

Finally, we consider the structural dynamics also in a fully time-dependent picture. Figure 6 presents the time-evolution of selected normal mode coordinates (i.e., those that display notable changes in amplitude over time and/or clear coherence) over the whole ensemble of trajectories. To gain some measure of the coherence between trajectories, we compare the mean values of the signed and unsigned coordinates. Equal amplitudes (of, e.g., a periodic oscillation or systematic
<figure><img src="image_10.png" /><figcaption>**FIG. 6.** Time evolution of selected normal mode coordinate projections $p_i$ (as described in Fig. 5). The blue colored heat plots represent the binned trajectory density for a given time. Orange lines mark the ensemble average and black lines the ensemble average of the absolute values. The normal mode coordinates are visualized to the right, as an orange distortion of the optimized geometry.</figcaption></figure>

drift) of the two then signal a completely coherent motion, whereas a vanishing average of the signed coordinates oppositely signals a complete dephasing.

Based on the described analysis, the following picture emerges: the dynamics are initially dominated by coherent in-plane H rocking, particularly in modes 17 and 22. As the $S_1$ state reaches its peak population and starts to decay, these modes dephase, as signaled by a rapid decrease in amplitude for the signed coordinates, whereas the amplitudes of the unsigned coordinates remain comparable. Concurrently, out-of-plane motion steadily increases over the first 100 fs and thereafter remains comparable throughout the remainder of the simulations. This is primarily seen from N-H centric modes 7 and 8 with unsigned values that drastically increase before flattening out. Due to the arbitrariness in the sign of the out-of-plane modes, however, the signed values appear to be artificially dephased. Within the same time frame, mode 30, which is essentially a pure N-H stretch, decreases notably in amplitude, i.e., the bond contracts from the $S_1$ state and on. Finally, an additional number of H rocking modes, in particular, 2, 20, and 21, display coherent motion over the whole time range. A partial dephasing, from the initially high coherence, occurs steadily over time, seen from a larger decrease in the sign than the unsigned coordinates, but with clear oscillations that remain even to the end of the simulations.

### D. Proton transfer

A complete ESPT is not found in any of the currently presented trajectories. Yet, large variations in the bonding structure are still observed, as demonstrated by the maximal $r_{N-H}$ value of 1.41 Å and the minimal $r_{S-H}$ value of 1.99 Å, as compared to their means (with standard deviations) of 1.03 (±0.06) Å and 2.78 (±0.22) Å. While limitations in the statistics still cannot be excluded as a possible cause, this suggests that the absence of ESPT, in the current work, might instead be the result of limitations in the computational model. Such limitations could include (but not necessarily be limited to) exclusion of solvent molecules (which might, in fact, mediate both the proton detachment and reattachment in the measurements), insufficient flexibility in the electronic structure method, or limitations imposed by the classical dynamics. Nevertheless, the observed electronic and structural dynamics still allow us to speculate about which aspects of the photo-dynamics could promote ESPT.

Starting from the Franck-Condon picture, it is unclear why any of the electronic excitations should promote ESPT, as they all appear to transfer electronic density from the sulfur to the nitrogen site, which should increase the basicity of the former and acidity of the latter. A Mulliken analysis, in fact, shows that the initial -0.46 net charge of the sulfur site increases notably to -0.04 ($S_1$), -0.25 ($S_2$), -0.19 ($T_1$), and -0.09 ($T_2$) upon excitation, with smaller effects for the nitrogen site. This clearly suggests that the ESPT must instead result from dynamical or external factors.

In terms of the structural dynamics that occur on short timescales (~100 fs), mode 22 is most likely to promote ESPT as it clearly rocks the N-H and C-S bonds with respect to each other, which could provide a channel for ultrafast intra-molecular ESPT before the dephasing of modes 17 and 22 sets in. Such a mechanism would be consistent with the weak signature of N-H cleavage previously observed using femtosecond RIXS.²⁰ It could, however, prove to be difficult to satisfactorily investigate using computational models, as it may require both a quantum description of the dynamics (to, for instance, account for
<figure><img src="image_11.png" /><figcaption>**FIG. 6.** Time evolution of selected normal mode coordinate projections $p_i$ (as described in Fig. 5). The blue colored heat plots represent the binned trajectory density for a given time. Orange lines mark the ensemble average and black lines the ensemble average of the absolute values. The normal mode coordinates are visualized to the right, as an orange distortion of the optimized geometry.</figcaption></figure>

drift) of the two then signal a completely coherent motion, whereas a vanishing average of the signed coordinates oppositely signals a complete dephasing.

Based on the described analysis, the following picture emerges: the dynamics are initially dominated by coherent in-plane H rocking, particularly in modes 17 and 22. As the $S_1$ state reaches its peak population and starts to decay, these modes dephase, as signaled by a rapid decrease in amplitude for the signed coordinates, whereas the amplitudes of the unsigned coordinates remain comparable. Concurrently, out-of-plane motion steadily increases over the first 100 fs and thereafter remains comparable throughout the remainder of the simulations. This is primarily seen from N-H centric modes 7 and 8 with unsigned values that drastically increase before flattening out. Due to the arbitrariness in the sign of the out-of-plane modes, however, the signed values appear to be artificially dephased. Within the same time frame, mode 30, which is essentially a pure N-H stretch, decreases notably in amplitude, i.e., the bond contracts from the $S_1$ state and on. Finally, an additional number of H rocking modes, in particular, 2, 20, and 21, display coherent motion over the whole time range. A partial dephasing, from the initially high coherence, occurs steadily over time, seen from a larger decrease in the sign than the unsigned coordinates, but with clear oscillations that remain even to the end of the simulations.

### D. Proton transfer

A complete ESPT is not found in any of the currently presented trajectories. Yet, large variations in the bonding structure are still observed, as demonstrated by the maximal $r_{N-H}$ value of 1.41 Å and the minimal $r_{S-H}$ value of 1.99 Å, as compared to their means (with standard deviations) of 1.03 (±0.06) Å and 2.78 (±0.22) Å. While limitations in the statistics still cannot be excluded as a possible cause, this suggests that the absence of ESPT, in the current work, might instead be the result of limitations in the computational model. Such limitations could include (but not necessarily be limited to) exclusion of solvent molecules (which might, in fact, mediate both the proton detachment and reattachment in the measurements), insufficient flexibility in the electronic structure method, or limitations imposed by the classical dynamics. Nevertheless, the observed electronic and structural dynamics still allow us to speculate about which aspects of the photo-dynamics could promote ESPT.

Starting from the Franck-Condon picture, it is unclear why any of the electronic excitations should promote ESPT, as they all appear to transfer electronic density from the sulfur to the nitrogen site, which should increase the basicity of the former and acidity of the latter. A Mulliken analysis, in fact, shows that the initial -0.46 net charge of the sulfur site increases notably to -0.04 ($S_1$), -0.25 ($S_2$), -0.19 ($T_1$), and -0.09 ($T_2$) upon excitation, with smaller effects for the nitrogen site. This clearly suggests that the ESPT must instead result from dynamical or external factors.

In terms of the structural dynamics that occur on short timescales (~100 fs), mode 22 is most likely to promote ESPT as it clearly rocks the N-H and C-S bonds with respect to each other, which could provide a channel for ultrafast intra-molecular ESPT before the dephasing of modes 17 and 22 sets in. Such a mechanism would be consistent with the weak signature of N-H cleavage previously observed using femtosecond RIXS.²⁰ It could, however, prove to be difficult to satisfactorily investigate using computational models, as it may require both a quantum description of the dynamics (to, for instance, account for
tunneling of the proton through a classically forbidden barrier) and
significant statistics to either confirm or reject. Modes 2, 20, and 21
could similarly participate in an intra-molecular transfer as they also
affect the S-H distance, but they all exhibit time periods of roughly
150 fs, i.e., comparable to the time point where out-of-plane motion
instead sets in to dominate.

In less than 100 fs, the N-H bond on average contracts, which
should inhibit its cleavage. ESPT on longer timescales is, therefore,
likely to be mediated by external factors, which strongly points toward
interactions with one or several solvent molecules. The bond contrac-
tion onsets together with the strong out-of-plane motions, in particu-
lar, N-H centered modes 7 and 8, which, in a solution environment,
could strongly affect the proton's coordination with the solvent. In our
previously mentioned Born-Oppenheimer molecular dynamics
study,²³ we further found both that the compound exhibits strong
hydrogen bonding with aqueous solutions and that the nitrogen site
coordinated water molecule of the anion species 2-TP⁻ interchanged
through a rotation that realigned the molecular plane with respect to
the solvent. A similar mechanism, i.e., solvent re-coordination through
out-of-plane motion, could, therefore, also be responsible for solvent-
mediated ESPT at timescales above 100 fs. This would be consistent
with pico- to nanosecond NEXAFS measurements, which find the
ESPT to be a rare event, presumably mediated by large fluctuations
that concurrently result in T₁ to S₀ relaxation.

Based on currently presented and recent insights,²³,²⁴ we finally
want to point out two clear ways that time-resolved spectroscopy
could be constructively employed to shed further light on the ESPT of
2-TP: (i) femtosecond timescale measurements, with a probe capable
of resolving either electronic states or vibrational modes, would be par-
ticularly useful to evaluate speculations about the ESPT mechanism
and pathway from current and previous work.¹⁸⁻²⁰ (ii) Gas phase mea-
surements on any timescale from femto- to pico- to nanoseconds,
with a probe sensitive to the nitrogen and/or sulfur protonation, could
conclusively determine whether a solvent is necessary to mediate the
ESPT. Insights derived from such experiments could, possibly, fill all
the significant gaps in the current understanding of the pathway to
fully “solve” the ESPT process of 2-TP.

IV. CONCLUSIONS

Non-adiabatic dynamics of 2-TP have been simulated by surface
hopping simulations including spin-orbit induced transitions, model-
ing the UV-induced excited-state dynamics probed in multiple previ-
ous time-resolved spectroscopic measurements.¹⁸,²⁰⁻²²

An ultrafast S₂ to S₁ population transfer, mediated by coherent
in-plane vibrations, is followed by complex population kinetics that
appreciably populates all the energetically available states through
multiple pathways. The complexity results from out-of-plane struc-
tural dynamics that onset in the S₁ state and significantly mixes the
excitation character of the electronic states, which, in particular, drives
ISC through the El-Sayed forbidden S₁(n, π⁺) to T₂(n, π⁺) channel.
The transfer of the population from singlets to triplets, once summed
over the individual states, can nevertheless be well described as a sim-
ple exponential decay with a lifetime of ~1 ps. This follows the trend
of efficient triplet population observed for multiple similar thio-
substituted compounds and offers rationalization for why the same
excited-stated dynamics have been previously observed for both S₂
and S₄ excitation of the system, on picosecond scales and longer.²²

While no ESPT was observed in the current work, two aspects of
the structural dynamics are suggested as possible mediators for it in
previously reported experiments: initial coherent in-plane motion for
ultrafast (~100 fs) intra-molecular transfer and out-of-plane motion
for solvent-mediated transfer at longer timescales. With this in mind,
we suggest two ways in which time-resolved spectroscopy could be
employed to address important remaining questions about the ESPT
pathway.

SUPPLEMENTARY MATERIAL

See the supplementary material for comparison to dynamics per-
formed with a larger basis set, estimation of state energy splittings in
the dynamics, comparison of electronic state populations in other state
representations, comparison to other kinetic model fits, visualization
of the MECP geometries, visualization of hop-inducing geometries for
backward transitions, visualization of normal modes of the electronic
ground state, and example input files for the OpenMolcas and
SHARC computations.

ACKNOWLEDGMENTS

This work was supported by the Swedish Research Council (Contract No. VR 2015-03956), the Carl Tryggers Foundation (Contract No. CTS18:285), and the Helmholtz Virtual Institute VI419 "Dynamic Pathways in Multidimensional Landscapes." The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) under allocation SNIC2019-1-9 at the Swedish National Supercomputer Center (NSC), the High Performance Computer Center North (HPC2N), and the Chalmers Centre for Computational Science and Engineering (C3SE). We thank Markus Kowalewski for helpful discussions and Sebastian Mai for assistance in using SHARC.

REFERENCES

<footnote n="1">¹R. Improta, F. Santoro, and L. Blancafort, "Quantum mechanical studies on the photophysics and the photochemistry of nucleic acids and nucleobases," *Chem. Rev.* **116**, 3540–3593 (2016).</footnote>

<footnote n="2">²C. Reichardt and C. E. Crespo-Hernández, “Ultrafast spin crossover in 4-thiothymidine in an ionic liquid,” *Chem. Commun.* **46**, 5963–5965 (2010).</footnote>

<footnote n="3">³C. Reichardt and C. E. Crespo-Hernández, "Room-temperature phosphorescence of the DNA monomer analogue 4-thiothymidine in aqueous solutions after UVA excitation," *J. Phys. Chem. Lett.* **1**, 2239–2243 (2010).</footnote>

<footnote n="4">⁴C. Reichardt, C. Guo, and C. E. Crespo-Hernández, "Excited-state dynamics in 6-thioguanosine from the femtosecond to microsecond time scale," *J. Phys. Chem. B* **115**, 3263–3270 (2011).</footnote>

<footnote n="5">⁵K. Taras-Goślińska, G. Burdziński, and G. Wenska, "Relaxation of the T1 excited state of 2-thiothymine, its riboside and deoxyriboside-enhanced nonradiative decay rate induced by sugar substituent," *J. Photochem. Photobiol. A* **275**, 89–95 (2014).</footnote>

<footnote n="6">⁶M. Pollum and C. E. Crespo-Hernández, "Communication: The dark singlet state as a doorway state in the ultrafast and efficient intersystem crossing dynamics in 2-thiothymine and 2-thiouracil," *J. Chem. Phys.* **140**, 071101 (2014).</footnote>

<footnote n="7">⁷M. Pollum, S. Jockusch, and C. E. Crespo-Hernández, "2,4-dithiothymine as a potent UVA chemotherapeutic agent," *J. Am. Chem. Soc.* **136**, 17930–17933 (2014).</footnote>

<footnote n="8">⁸M. Pollum, S. Jockusch, and C. E. Crespo-Hernández, "Increase in the photo-reactivity of uracil derivatives by doubling thionation," *Phys. Chem. Chem. Phys.* **17**, 27851–27861 (2015).</footnote>
tunneling of the proton through a classically forbidden barrier) and
significant statistics to either confirm or reject. Modes 2, 20, and 21
could similarly participate in an intra-molecular transfer as they also
affect the S-H distance, but they all exhibit time periods of roughly
150 fs, i.e., comparable to the time point where out-of-plane motion
instead sets in to dominate.

In less than 100 fs, the N-H bond on average contracts, which
should inhibit its cleavage. ESPT on longer timescales is, therefore,
likely to be mediated by external factors, which strongly points toward
interactions with one or several solvent molecules. The bond contrac-
tion onsets together with the strong out-of-plane motions, in particu-
lar, N-H centered modes 7 and 8, which, in a solution environment,
could strongly affect the proton's coordination with the solvent. In our
previously mentioned Born-Oppenheimer molecular dynamics
study,²³ we further found both that the compound exhibits strong
hydrogen bonding with aqueous solutions and that the nitrogen site
coordinated water molecule of the anion species 2-TP⁻ interchanged
through a rotation that realigned the molecular plane with respect to
the solvent. A similar mechanism, i.e., solvent re-coordination through
out-of-plane motion, could, therefore, also be responsible for solvent-
mediated ESPT at timescales above 100 fs. This would be consistent
with pico- to nanosecond NEXAFS measurements, which find the
ESPT to be a rare event, presumably mediated by large fluctuations
that concurrently result in T₁ to S₀ relaxation.

Based on currently presented and recent insights,²³,²⁴ we finally
want to point out two clear ways that time-resolved spectroscopy
could be constructively employed to shed further light on the ESPT of
2-TP: (i) femtosecond timescale measurements, with a probe capable
of resolving either electronic states or vibrational modes, would be par-
ticularly useful to evaluate speculations about the ESPT mechanism
and pathway from current and previous work.¹⁸⁻²⁰ (ii) Gas phase mea-
surements on any timescale from femto- to pico- to nanoseconds,
with a probe sensitive to the nitrogen and/or sulfur protonation, could
conclusively determine whether a solvent is necessary to mediate the
ESPT. Insights derived from such experiments could, possibly, fill all
the significant gaps in the current understanding of the pathway to
fully “solve” the ESPT process of 2-TP.

IV. CONCLUSIONS

Non-adiabatic dynamics of 2-TP have been simulated by surface
hopping simulations including spin-orbit induced transitions, model-
ing the UV-induced excited-state dynamics probed in multiple previ-
ous time-resolved spectroscopic measurements.¹⁸,²⁰⁻²²

An ultrafast S₂ to S₁ population transfer, mediated by coherent
in-plane vibrations, is followed by complex population kinetics that
appreciably populates all the energetically available states through
multiple pathways. The complexity results from out-of-plane struc-
tural dynamics that onset in the S₁ state and significantly mixes the
excitation character of the electronic states, which, in particular, drives
ISC through the El-Sayed forbidden S₁(n, π⁺) to T₂(n, π⁺) channel.
The transfer of the population from singlets to triplets, once summed
over the individual states, can nevertheless be well described as a sim-
ple exponential decay with a lifetime of ~1 ps. This follows the trend
of efficient triplet population observed for multiple similar thio-
substituted compounds and offers rationalization for why the same
excited-stated dynamics have been previously observed for both S₂
and S₄ excitation of the system, on picosecond scales and longer.²²

While no ESPT was observed in the current work, two aspects of
the structural dynamics are suggested as possible mediators for it in
previously reported experiments: initial coherent in-plane motion for
ultrafast (~100 fs) intra-molecular transfer and out-of-plane motion
for solvent-mediated transfer at longer timescales. With this in mind,
we suggest two ways in which time-resolved spectroscopy could be
employed to address important remaining questions about the ESPT
pathway.

SUPPLEMENTARY MATERIAL

See the supplementary material for comparison to dynamics per-
formed with a larger basis set, estimation of state energy splittings in
the dynamics, comparison of electronic state populations in other state
representations, comparison to other kinetic model fits, visualization
of the MECP geometries, visualization of hop-inducing geometries for
backward transitions, visualization of normal modes of the electronic
ground state, and example input files for the OpenMolcas and
SHARC computations.

ACKNOWLEDGMENTS

This work was supported by the Swedish Research Council (Contract No. VR 2015-03956), the Carl Tryggers Foundation (Contract No. CTS18:285), and the Helmholtz Virtual Institute VI419 "Dynamic Pathways in Multidimensional Landscapes." The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) under allocation SNIC2019-1-9 at the Swedish National Supercomputer Center (NSC), the High Performance Computer Center North (HPC2N), and the Chalmers Centre for Computational Science and Engineering (C3SE). We thank Markus Kowalewski for helpful discussions and Sebastian Mai for assistance in using SHARC.

REFERENCES

<footnote n="1">¹R. Improta, F. Santoro, and L. Blancafort, "Quantum mechanical studies on the photophysics and the photochemistry of nucleic acids and nucleobases," *Chem. Rev.* **116**, 3540–3593 (2016).</footnote>

<footnote n="2">²C. Reichardt and C. E. Crespo-Hernández, “Ultrafast spin crossover in 4-thiothymidine in an ionic liquid,” *Chem. Commun.* **46**, 5963–5965 (2010).</footnote>

<footnote n="3">³C. Reichardt and C. E. Crespo-Hernández, "Room-temperature phosphorescence of the DNA monomer analogue 4-thiothymidine in aqueous solutions after UVA excitation," *J. Phys. Chem. Lett.* **1**, 2239–2243 (2010).</footnote>

<footnote n="4">⁴C. Reichardt, C. Guo, and C. E. Crespo-Hernández, "Excited-state dynamics in 6-thioguanosine from the femtosecond to microsecond time scale," *J. Phys. Chem. B* **115**, 3263–3270 (2011).</footnote>

<footnote n="5">⁵K. Taras-Goślińska, G. Burdziński, and G. Wenska, "Relaxation of the T1 excited state of 2-thiothymine, its riboside and deoxyriboside-enhanced nonradiative decay rate induced by sugar substituent," *J. Photochem. Photobiol. A* **275**, 89–95 (2014).</footnote>

<footnote n="6">⁶M. Pollum and C. E. Crespo-Hernández, "Communication: The dark singlet state as a doorway state in the ultrafast and efficient intersystem crossing dynamics in 2-thiothymine and 2-thiouracil," *J. Chem. Phys.* **140**, 071101 (2014).</footnote>

<footnote n="7">⁷M. Pollum, S. Jockusch, and C. E. Crespo-Hernández, "2,4-dithiothymine as a potent UVA chemotherapeutic agent," *J. Am. Chem. Soc.* **136**, 17930–17933 (2014).</footnote>

<footnote n="8">⁸M. Pollum, S. Jockusch, and C. E. Crespo-Hernández, "Increase in the photo-reactivity of uracil derivatives by doubling thionation," *Phys. Chem. Chem. Phys.* **17**, 27851–27861 (2015).</footnote>
<footnote n="9">⁹M. Pollum, L. A. Ortiz-Rodríguez, S. Jockusch, and C. E. Crespo-Hernández, “The triplet state of 6-thio-2-deoxyguanosine: Intrinsic properties and reactivity toward molecular oxygen,” *Photochem. Photobiol.* **92**, 286–292 (2016).</footnote>

<footnote n="10">¹⁰L. Martínez-Fernández, L. González, and I. Corral, “An ab initio mechanism for efficient population of triplet states in cytotoxic sulfur substituted DNA bases: The case of 6-thioguanine,” *Chem. Commun.* **48**, 2134–2136 (2012).</footnote>

<footnote n="11">¹¹G. Cui and W.-H. Fang, “State-specific heavy-atom effect on intersystem crossing processes in 2-thiothymine: A potential photodynamic therapy photosensitizer,” *J. Chem. Phys.* **138**, 044315 (2013).</footnote>

<footnote n="12">¹²J. P. Gobbo and A. C. Borin, “2-Thiouracil deactivation pathways and triplet states population,” *Comput. Theor. Chem.* **1040-1041**, 195–201 (2014).</footnote>

<footnote n="13">¹³G. Cui and W. Thiel, “Intersystem crossing enables 4-thiothymidine to act as a photosensitizer in photodynamic therapy: An ab initio QM/MM study,” *J. Phys. Chem. Lett.* **5**, 2682–2687 (2014).</footnote>

<footnote n="14">¹⁴L. Martínez-Fernández, I. Corral, G. Granucci, and M. Persico, “Competing ultrafast intersystem crossing and internal conversion: A time resolved picture for the deactivation of 6-thioguanine,” *Chem. Sci.* **5**, 1336–1347 (2014).</footnote>

<footnote n="15">¹⁵J. Pirillo, B. C. De Simone, and N. Russo, “Photophysical properties prediction of selenium- and tellurium-substituted thymidine as potential UVA chemo-therapeutic agents,” *Theor. Chem. Acc.* **135**, 8 (2016).</footnote>

<footnote n="16">¹⁶H. Yu, J. A. Sanchez-Rodriguez, M. Pollum, C. E. Crespo-Hernández, S. Mai, P. Marquetand, L. González, and S. Ullrich, “Internal conversion and intersystem crossing pathways in UV excited, isolated uracils and their implications in prebiotic chemistry,” *Phys. Chem. Chem. Phys.* **18**, 20168–20176 (2016).</footnote>

<footnote n="17">¹⁷S. Mai, P. Marquetand, and L. González, “Intersystem crossing pathways in the noncanonical nucleobase 2-thiouracil: A time-dependent picture,” *J. Phys. Chem. Lett.* **7**, 1978–1983 (2016).</footnote>

<footnote n="18">¹⁸R. Du, C. Liu, Y. Zhao, K.-M. Pei, H.-G. Wang, X. Zheng, M. Li, J.-D. Xue, and D. L. Phillips, “Resonance Raman spectroscopic and theoretical investigation of the excited state proton transfer reaction dynamics of 2-thiopyridone,” *J. Phys. Chem. B* **115**, 8266–8277 (2011).</footnote>

<footnote n="19">¹⁹S. Eckert, P. Miedema, W. Quevedo, B. O’Cinneide, M. Fondell, M. Beye, A. Pietzsch, M. Ross, M. Khalil, and A. Föhlisch, “Molecular structures and protonation state of 2-Mercaptopyridine in aqueous solution,” *Chem. Phys. Lett.* **647**, 103–106 (2016).</footnote>

<footnote n="20">²⁰S. Eckert, J. Norell, P. S. Miedema, M. Beye, M. Fondell, W. Quevedo, B. Kennedy, M. Hantschmann, A. Pietzsch, B. E. Van Kuiken, M. Ross, M. P. Minitti, S. P. Moeller, W. F. Schlotter, M. Khalil, M. Odelius, and A. Föhlisch, “Ultrafast independent N-H and N-C bond deformation investigated with resonant inelastic x-ray scattering,” *Angew. Chem., Int. Ed.* **56**, 6088–6092 (2017).</footnote>

<footnote n="21">²¹B. E. Van Kuiken, M. R. Ross, M. L. Strader, A. A. Cordones, H. Cho, J. H. Lee, R. W. Schoenlein, and M. Khalil, “Picosecond sulfur K-edge X-ray absorption spectroscopy with applications to excited state proton transfer,” *Struct. Dyn.* **4**, 044021 (2017).</footnote>

<footnote n="22">²²S. Eckert, J. Norell, R. M. Jay, M. Fondell, R. Mitzner, M. Odelius, and A. Föhlisch, “T1 population as the driver of excited-state proton-transfer in 2-thiopyridone,” *Chemistry* **25**, 1733–1739 (2019).</footnote>

<footnote n="23">²³J. Norell, A. Ljungdahl, and M. Odelius, “Interdependent electronic structure, protonation, and solvatization of aqueous 2-thiopyridone,” *J. Phys. Chem. B* **123**, 5555–5567 (2019).</footnote>

<footnote n="24">²⁴J. Norell, S. Eckert, B. E. Van Kuiken, A. Föhlisch, and M. Odelius, “Ab initio simulations of complementary k-edges and solvatization effects for detection of proton transfer in aqueous 2-thiopyridone,” *J. Chem. Phys.* **151**, 114117 (2019).</footnote>

<footnote n="25">²⁵Y. Zhang, K. de La Harpe, A. A. Beckstead, R. Improta, and B. Kohler, “UV-induced proton transfer between DNA strands,” *J. Am. Chem. Soc.* **137**, 7059–7062 (2015).</footnote>

<footnote n="26">²⁶A. Corani, A. Huijser, T. Gustavsson, D. Markovitsi, P.-Å. Malmqvist, A. Pezzella, M. d’Ischia, and V. Sundström, “Superior photoprotective motifs and mechanisms in eumelanins uncovered,” *J. Am. Chem. Soc.* **136**, 11626–11635 (2014).</footnote>

<footnote n="27">²⁷R. Jones and A. Katritzky, “721. Tautomeric pyridines. Part I. Pyrid-2-and-4-thione,” *J. Chem. Soc.* **1958**, 3610–3613 (1958).</footnote>

<footnote n="28">²⁸E. Spinner, “250. The infrared spectra of some N-heteroaromatic mercapto-compounds and of their N-methyl and S-methyl derivatives,” *J. Chem. Soc.* **1960**, 1237–1242 (1960).</footnote>

<footnote n="29">²⁹P. Beak, F. S. Fry, Jr., J. Lee, and F. Steele, “Equilibration studies. Protomeric equilibria of 2- and 4-hydroxypyridines, 2- and 4-hydroxypyrimidines, 2- and 4-mercaptopyridines, and structurally related compounds in the gas phase,” *J. Am. Chem. Soc.* **98**, 171–179 (1976).</footnote>

<footnote n="30">³⁰S. Stoyanov, I. Petkov, L. Antonov, T. Stoyanova, P. Karagiannidis, and P. Aslanidis, “Thione-thiol tautomerism and stability of 2- and 4-mercaptopyridines and 2-mercaptopyrimidines,” *Can. J. Chem.* **68**, 1482–1489 (1990).</footnote>

<footnote n="31">³¹S. Melandri, L. Evangelisti, A. Maris, W. Caminati, B. M. Giuliano, V. Feyer, K. C. Prince, and M. Coreno, “Rotational and core level spectroscopies as complementary techniques in tautomeric/conformational studies: The case of 2-mercaptopyridine,” *J. Am. Chem. Soc.* **132**, 10269–10271 (2010).</footnote>

<footnote n="32">³²S. Mai, M. Pollum, L. Martínez-Fernández, N. Dunn, P. Marquetand, I. Corral, C. E. Crespo-Hernández, and L. González, “The origin of efficient triplet state population in sulfur-substituted nucleobases,” *Nat. Commun.* **7**, 13077 (2016).</footnote>

<footnote n="33">³³J. C. Tully, “Molecular dynamics with electronic transitions,” *J. Chem. Phys.* **93**, 1061–1071 (1990).</footnote>

<footnote n="34">³⁴M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, and L. González, “Sharc: Ab initio molecular dynamics with surface hopping in the adiabatic representation including arbitrary couplings,” *J. Chem. Theory Comput.* **7**, 1253–1258 (2011).</footnote>

<footnote n="35">³⁵S. Mai, P. Marquetand, and L. González, “Nonadiabatic dynamics: The SHARC approach,” *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **8**, e1370 (2018).</footnote>

<footnote n="36">³⁶S. Mai, P. Marquetand, and L. González, “A general method to describe intersystem crossing dynamics in trajectory surface hopping,” *Int. J. Quantum Chem.* **115**, 1215–1231 (2015).</footnote>

<footnote n="37">³⁷G. Granucci, M. Persico, and A. Toniolo, “Direct semiclassical simulation of photochemical processes with semiempirical wave functions,” *J. Chem. Phys.* **114**, 10608–10615 (2001).</footnote>

<footnote n="38">³⁸J. Pittner, H. Lischka, and M. Barbatti, “Optimization of mixed quantum-classical dynamics: Time-derivative coupling terms and selected couplings,” *Chem. Phys.* **356**, 147–152 (2009).</footnote>

<footnote n="39">³⁹F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico, and H. Lischka, “Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer,” *J. Chem. Phys.* **137**, 22A514 (2012).</footnote>

<footnote n="40">⁴⁰G. Granucci and M. Persico, “Critical appraisal of the fewest switches algorithm for surface hopping,” *J. Chem. Phys.* **126**, 134114 (2007).</footnote>

<footnote n="41">⁴¹I. Fdez. Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov, R. K. Carlson, L. F. Chibotaru, J. Creutzberg, N. Dattani, M. G. Delcey, S. S. Dong, A. Dreuw, L. Freitag, L. M. Frutos, L. Gagliardi, F. Gendron, A. Giussani, L. González, G. Grell, M. Guo, C. E. Hoyer, M. Johansson, S. Keller, S. Knecht, G. Kovacevic, E. Källman, G. L. Manni, M. Lundberg, Y. Ma, S. Mai, J. P. Malhado, P. Å. Malmqvist, P. Marquetand, S. A. Mewes, J. Norell, M. Olivucci, M. Oppel, Q. M. Phung, K. Pierloot, F. Plasser, M. Reiher, A. M. Sand, I. Schapiro, P. Sharma, C. J. Stein, L. K. Sørensen, D. G. Truhlar, M. Ugandi, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, O. Weser, T. A. Wesolowski, P.-O. Widmark, S. Wouters, A. Zech, J. P. Zobel, and R. Lindh, “OpenMolcas: From source code to insight,” *J. Chem. Theory Comput.* **15**, 5925–5964 (2019).</footnote>

<footnote n="42">⁴²B. O. Roos, P. R. Taylor, and P. E. Siegbahn, “A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach,” *Chem. Phys.* **48**, 157–173 (1980).</footnote>

<footnote n="43">⁴³K. Andersson, P.-Å. Malmqvist, and B. O. Roos, “Second-order perturbation theory with a complete active space self-consistent field reference function,” *J. Chem. Phys.* **96**, 1218–1226 (1992).</footnote>

<footnote n="44">⁴⁴J. Finley, P.-Å. Malmqvist, B. O. Roos, and L. Serrano-Andrés, “The multi-state CASPT2 method,” *Chem. Phys. Lett.* **288**, 299–306 (1998).</footnote>

<footnote n="45">⁴⁵M. Douglas and N. M. Kroll, “Quantum electrodynamical corrections to the fine structure of helium,” *Ann. Phys.* **82**, 89–155 (1974).</footnote>

<footnote n="46">⁴⁶B. A. Hess, “Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators,” *Phys. Rev. A* **33**, 3742–3748 (1986).</footnote>

<footnote n="47">⁴⁷B. O. Roos, R. Lindh, P. Å. Malmqvist, V. Veryazov, and P. O. Widmark, “Main group atoms and dimers studied with a new relativistic ANO basis set,” *J. Phys. Chem. A* **108**, 2851–2858 (2004).</footnote>
<footnote n="9">⁹M. Pollum, L. A. Ortiz-Rodríguez, S. Jockusch, and C. E. Crespo-Hernández, “The triplet state of 6-thio-2-deoxyguanosine: Intrinsic properties and reactivity toward molecular oxygen,” *Photochem. Photobiol.* **92**, 286–292 (2016).</footnote>

<footnote n="10">¹⁰L. Martínez-Fernández, L. González, and I. Corral, “An ab initio mechanism for efficient population of triplet states in cytotoxic sulfur substituted DNA bases: The case of 6-thioguanine,” *Chem. Commun.* **48**, 2134–2136 (2012).</footnote>

<footnote n="11">¹¹G. Cui and W.-H. Fang, “State-specific heavy-atom effect on intersystem crossing processes in 2-thiothymine: A potential photodynamic therapy photosensitizer,” *J. Chem. Phys.* **138**, 044315 (2013).</footnote>

<footnote n="12">¹²J. P. Gobbo and A. C. Borin, “2-Thiouracil deactivation pathways and triplet states population,” *Comput. Theor. Chem.* **1040-1041**, 195–201 (2014).</footnote>

<footnote n="13">¹³G. Cui and W. Thiel, “Intersystem crossing enables 4-thiothymidine to act as a photosensitizer in photodynamic therapy: An ab initio QM/MM study,” *J. Phys. Chem. Lett.* **5**, 2682–2687 (2014).</footnote>

<footnote n="14">¹⁴L. Martínez-Fernández, I. Corral, G. Granucci, and M. Persico, “Competing ultrafast intersystem crossing and internal conversion: A time resolved picture for the deactivation of 6-thioguanine,” *Chem. Sci.* **5**, 1336–1347 (2014).</footnote>

<footnote n="15">¹⁵J. Pirillo, B. C. De Simone, and N. Russo, “Photophysical properties prediction of selenium- and tellurium-substituted thymidine as potential UVA chemo-therapeutic agents,” *Theor. Chem. Acc.* **135**, 8 (2016).</footnote>

<footnote n="16">¹⁶H. Yu, J. A. Sanchez-Rodriguez, M. Pollum, C. E. Crespo-Hernández, S. Mai, P. Marquetand, L. González, and S. Ullrich, “Internal conversion and intersystem crossing pathways in UV excited, isolated uracils and their implications in prebiotic chemistry,” *Phys. Chem. Chem. Phys.* **18**, 20168–20176 (2016).</footnote>

<footnote n="17">¹⁷S. Mai, P. Marquetand, and L. González, “Intersystem crossing pathways in the noncanonical nucleobase 2-thiouracil: A time-dependent picture,” *J. Phys. Chem. Lett.* **7**, 1978–1983 (2016).</footnote>

<footnote n="18">¹⁸R. Du, C. Liu, Y. Zhao, K.-M. Pei, H.-G. Wang, X. Zheng, M. Li, J.-D. Xue, and D. L. Phillips, “Resonance Raman spectroscopic and theoretical investigation of the excited state proton transfer reaction dynamics of 2-thiopyridone,” *J. Phys. Chem. B* **115**, 8266–8277 (2011).</footnote>

<footnote n="19">¹⁹S. Eckert, P. Miedema, W. Quevedo, B. O’Cinneide, M. Fondell, M. Beye, A. Pietzsch, M. Ross, M. Khalil, and A. Föhlisch, “Molecular structures and protonation state of 2-Mercaptopyridine in aqueous solution,” *Chem. Phys. Lett.* **647**, 103–106 (2016).</footnote>

<footnote n="20">²⁰S. Eckert, J. Norell, P. S. Miedema, M. Beye, M. Fondell, W. Quevedo, B. Kennedy, M. Hantschmann, A. Pietzsch, B. E. Van Kuiken, M. Ross, M. P. Minitti, S. P. Moeller, W. F. Schlotter, M. Khalil, M. Odelius, and A. Föhlisch, “Ultrafast independent N-H and N-C bond deformation investigated with resonant inelastic x-ray scattering,” *Angew. Chem., Int. Ed.* **56**, 6088–6092 (2017).</footnote>

<footnote n="21">²¹B. E. Van Kuiken, M. R. Ross, M. L. Strader, A. A. Cordones, H. Cho, J. H. Lee, R. W. Schoenlein, and M. Khalil, “Picosecond sulfur K-edge X-ray absorption spectroscopy with applications to excited state proton transfer,” *Struct. Dyn.* **4**, 044021 (2017).</footnote>

<footnote n="22">²²S. Eckert, J. Norell, R. M. Jay, M. Fondell, R. Mitzner, M. Odelius, and A. Föhlisch, “T1 population as the driver of excited-state proton-transfer in 2-thiopyridone,” *Chemistry* **25**, 1733–1739 (2019).</footnote>

<footnote n="23">²³J. Norell, A. Ljungdahl, and M. Odelius, “Interdependent electronic structure, protonation, and solvatization of aqueous 2-thiopyridone,” *J. Phys. Chem. B* **123**, 5555–5567 (2019).</footnote>

<footnote n="24">²⁴J. Norell, S. Eckert, B. E. Van Kuiken, A. Föhlisch, and M. Odelius, “Ab initio simulations of complementary k-edges and solvatization effects for detection of proton transfer in aqueous 2-thiopyridone,” *J. Chem. Phys.* **151**, 114117 (2019).</footnote>

<footnote n="25">²⁵Y. Zhang, K. de La Harpe, A. A. Beckstead, R. Improta, and B. Kohler, “UV-induced proton transfer between DNA strands,” *J. Am. Chem. Soc.* **137**, 7059–7062 (2015).</footnote>

<footnote n="26">²⁶A. Corani, A. Huijser, T. Gustavsson, D. Markovitsi, P.-Å. Malmqvist, A. Pezzella, M. d’Ischia, and V. Sundström, “Superior photoprotective motifs and mechanisms in eumelanins uncovered,” *J. Am. Chem. Soc.* **136**, 11626–11635 (2014).</footnote>

<footnote n="27">²⁷R. Jones and A. Katritzky, “721. Tautomeric pyridines. Part I. Pyrid-2-and-4-thione,” *J. Chem. Soc.* **1958**, 3610–3613 (1958).</footnote>

<footnote n="28">²⁸E. Spinner, “250. The infrared spectra of some N-heteroaromatic mercapto-compounds and of their N-methyl and S-methyl derivatives,” *J. Chem. Soc.* **1960**, 1237–1242 (1960).</footnote>

<footnote n="29">²⁹P. Beak, F. S. Fry, Jr., J. Lee, and F. Steele, “Equilibration studies. Protomeric equilibria of 2- and 4-hydroxypyridines, 2- and 4-hydroxypyrimidines, 2- and 4-mercaptopyridines, and structurally related compounds in the gas phase,” *J. Am. Chem. Soc.* **98**, 171–179 (1976).</footnote>

<footnote n="30">³⁰S. Stoyanov, I. Petkov, L. Antonov, T. Stoyanova, P. Karagiannidis, and P. Aslanidis, “Thione-thiol tautomerism and stability of 2- and 4-mercaptopyridines and 2-mercaptopyrimidines,” *Can. J. Chem.* **68**, 1482–1489 (1990).</footnote>

<footnote n="31">³¹S. Melandri, L. Evangelisti, A. Maris, W. Caminati, B. M. Giuliano, V. Feyer, K. C. Prince, and M. Coreno, “Rotational and core level spectroscopies as complementary techniques in tautomeric/conformational studies: The case of 2-mercaptopyridine,” *J. Am. Chem. Soc.* **132**, 10269–10271 (2010).</footnote>

<footnote n="32">³²S. Mai, M. Pollum, L. Martínez-Fernández, N. Dunn, P. Marquetand, I. Corral, C. E. Crespo-Hernández, and L. González, “The origin of efficient triplet state population in sulfur-substituted nucleobases,” *Nat. Commun.* **7**, 13077 (2016).</footnote>

<footnote n="33">³³J. C. Tully, “Molecular dynamics with electronic transitions,” *J. Chem. Phys.* **93**, 1061–1071 (1990).</footnote>

<footnote n="34">³⁴M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, and L. González, “Sharc: Ab initio molecular dynamics with surface hopping in the adiabatic representation including arbitrary couplings,” *J. Chem. Theory Comput.* **7**, 1253–1258 (2011).</footnote>

<footnote n="35">³⁵S. Mai, P. Marquetand, and L. González, “Nonadiabatic dynamics: The SHARC approach,” *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **8**, e1370 (2018).</footnote>

<footnote n="36">³⁶S. Mai, P. Marquetand, and L. González, “A general method to describe intersystem crossing dynamics in trajectory surface hopping,” *Int. J. Quantum Chem.* **115**, 1215–1231 (2015).</footnote>

<footnote n="37">³⁷G. Granucci, M. Persico, and A. Toniolo, “Direct semiclassical simulation of photochemical processes with semiempirical wave functions,” *J. Chem. Phys.* **114**, 10608–10615 (2001).</footnote>

<footnote n="38">³⁸J. Pittner, H. Lischka, and M. Barbatti, “Optimization of mixed quantum-classical dynamics: Time-derivative coupling terms and selected couplings,” *Chem. Phys.* **356**, 147–152 (2009).</footnote>

<footnote n="39">³⁹F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico, and H. Lischka, “Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer,” *J. Chem. Phys.* **137**, 22A514 (2012).</footnote>

<footnote n="40">⁴⁰G. Granucci and M. Persico, “Critical appraisal of the fewest switches algorithm for surface hopping,” *J. Chem. Phys.* **126**, 134114 (2007).</footnote>

<footnote n="41">⁴¹I. Fdez. Galván, M. Vacher, A. Alavi, C. Angeli, F. Aquilante, J. Autschbach, J. J. Bao, S. I. Bokarev, N. A. Bogdanov, R. K. Carlson, L. F. Chibotaru, J. Creutzberg, N. Dattani, M. G. Delcey, S. S. Dong, A. Dreuw, L. Freitag, L. M. Frutos, L. Gagliardi, F. Gendron, A. Giussani, L. González, G. Grell, M. Guo, C. E. Hoyer, M. Johansson, S. Keller, S. Knecht, G. Kovacevic, E. Källman, G. L. Manni, M. Lundberg, Y. Ma, S. Mai, J. P. Malhado, P. Å. Malmqvist, P. Marquetand, S. A. Mewes, J. Norell, M. Olivucci, M. Oppel, Q. M. Phung, K. Pierloot, F. Plasser, M. Reiher, A. M. Sand, I. Schapiro, P. Sharma, C. J. Stein, L. K. Sørensen, D. G. Truhlar, M. Ugandi, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, O. Weser, T. A. Wesolowski, P.-O. Widmark, S. Wouters, A. Zech, J. P. Zobel, and R. Lindh, “OpenMolcas: From source code to insight,” *J. Chem. Theory Comput.* **15**, 5925–5964 (2019).</footnote>

<footnote n="42">⁴²B. O. Roos, P. R. Taylor, and P. E. Siegbahn, “A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach,” *Chem. Phys.* **48**, 157–173 (1980).</footnote>

<footnote n="43">⁴³K. Andersson, P.-Å. Malmqvist, and B. O. Roos, “Second-order perturbation theory with a complete active space self-consistent field reference function,” *J. Chem. Phys.* **96**, 1218–1226 (1992).</footnote>

<footnote n="44">⁴⁴J. Finley, P.-Å. Malmqvist, B. O. Roos, and L. Serrano-Andrés, “The multi-state CASPT2 method,” *Chem. Phys. Lett.* **288**, 299–306 (1998).</footnote>

<footnote n="45">⁴⁵M. Douglas and N. M. Kroll, “Quantum electrodynamical corrections to the fine structure of helium,” *Ann. Phys.* **82**, 89–155 (1974).</footnote>

<footnote n="46">⁴⁶B. A. Hess, “Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators,” *Phys. Rev. A* **33**, 3742–3748 (1986).</footnote>

<footnote n="47">⁴⁷B. O. Roos, R. Lindh, P. Å. Malmqvist, V. Veryazov, and P. O. Widmark, “Main group atoms and dimers studied with a new relativistic ANO basis set,” *J. Phys. Chem. A* **108**, 2851–2858 (2004).</footnote>
<footnote n="48">⁴⁸F. Aquilante, T. B. Pedersen, and R. Lindh, "Low-cost evaluation of the exchange Fock matrix from Cholesky and density fitting representations of the electron repulsion integrals," *J. Chem. Phys.* **126**, 194106 (2007).</footnote>

<footnote n="49">⁴⁹J. Boström, M. G. Delcey, F. Aquilante, L. Serrano-Andrés, T. B. Pedersen, and R. Lindh, "Calibration of Cholesky auxiliary basis sets for multiconfigurational perturbation theory calculations of excitation energies," *J. Chem. Theory Comput.* **6**, 747–754 (2010).</footnote>

<footnote n="50">⁵⁰B. A. Hess, C. M. Marian, U. Wahlgren, and O. Gropen, "A mean-field spin-orbit method applicable to correlated wavefunctions," *Chem. Phys. Lett.* **251**, 365–371 (1996).</footnote>

<footnote n="51">⁵¹P.-À. Malmqvist and B. O. Roos, "The CASSCF state interaction method," *Chem. Phys. Lett.* **155**, 189–194 (1989).</footnote>

<footnote n="52">⁵²P.-À. Malmqvist, B. O. Roos, and B. Schimmelpfennig, "The restricted active space (RAS) state interaction approach with spin-orbit coupling," *Chem. Phys. Lett.* **357**, 230–240 (2002).</footnote>

<footnote n="53">⁵³G. Kovacevic and V. Veryazov, "Luscus: Molecular viewer and editor for MOLCAS," *J. Cheminf.* **7**, 16 (2015).</footnote>

<footnote n="54">⁵⁴T. J. Martínez, M. Ben-Nun, and R. D. Levine, "Molecular collision dynamics on several electronic states," *J. Phys. Chem. A* **101**, 6389–6402 (1997).</footnote>

<footnote n="55">⁵⁵B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martínez, "Implementation of ab initio multiple spawning in the Molpro quantum chemistry package," *Chem. Phys.* **347**, 3–16 (2008).</footnote>

<footnote n="56">⁵⁶J. D. Coe, B. G. Levine, and T. J. Martínez, "Ab initio molecular dynamics of excited-state intramolecular proton transfer using multireference perturbation theory," *J. Phys. Chem. A* **111**, 11302–11310 (2007).</footnote>

<footnote n="57">⁵⁷H. Tao, B. G. Levine, and T. J. Martínez, "Ab initio multiple spawning dynamics using multi-state second-order perturbation theory," *J. Phys. Chem. A* **113**, 13656–13662 (2009).</footnote>

<footnote n="58">⁵⁸H. Tao, T. K. Allison, T. W. Wright, A. M. Stooke, C. Khurmi, J. van Tilborg, Y. Liu, R. W. Falcone, A. Belkacem, and T. J. Martinez, "Ultrafast internal conversion in ethylene. I. The excited state lifetime," *J. Chem. Phys.* **134**, 244306 (2011).</footnote>

<footnote n="59">⁵⁹T. Mori, W. J. Glover, M. S. Schuurman, and T. J. Martinez, "Role of Rydberg states in the photochemical dynamics of ethylene," *J. Phys. Chem. A* **116**, 2808–2818 (2012).</footnote>

<footnote n="60">⁶⁰T. S. Kuhlman, W. J. Glover, T. Mori, K. B. Møller, and T. J. Martínez, "Between ethylene and polyenes—the non-adiabatic dynamics of cis-dienes," *Faraday Discuss.* **157**, 193–212 (2012).</footnote>

<footnote n="61">⁶¹T. J. A. Wolf, T. S. Kuhlman, O. Schalk, T. J. Martínez, K. B. Møller, A. Stolow, and A.-N. Unterreiner, "Hexamethylcyclopentadiene: Time-resolved photoelectron spectroscopy and ab initio multiple spawning simulations," *Phys. Chem. Chem. Phys.* **16**, 11770–11779 (2014).</footnote>

<footnote n="62">⁶²L. Liu, J. Liu, and T. J. Martínez, "Dynamical correlation effects on photoisomerization: Ab initio multiple spawning dynamics with MS-CASPT2 for a model trans-protonated Schiff base," *J. Phys. Chem. B* **120**, 1940–1949 (2016).</footnote>

<footnote n="63">⁶³W. J. Glover, T. Mori, M. S. Schuurman, A. E. Boguslavskiy, O. Schalk, A. Stolow, and T. J. Martínez, "Excited state non-adiabatic dynamics of the smallest polyene, trans 1,3-butadiene. II. Ab initio multiple spawning simulations," *J. Chem. Phys.* **148**, 164303 (2018).</footnote>

<footnote n="64">⁶⁴J. W. Park and T. Shiozaki, "On-the-fly CASPT2 surface-hopping dynamics," *J. Chem. Theory Comput.* **13**, 3676–3683 (2017).</footnote>

<footnote n="65">⁶⁵M. Richter, S. Mai, P. Marquetand, and L. González, "Ultrafast intersystem crossing dynamics in uracil unravelled by ab initio molecular dynamics," *Phys. Chem. Chem. Phys.* **16**, 24423–24436 (2014).</footnote>

<footnote n="66">⁶⁶S. Mai, P. Marquetand, M. Richter, J. González-Vázquez, and L. González, "Singlet and triplet excited-state dynamics study of the keto and enol tautomers of cytosine," *ChemPhysChem* **14**, 2920–2931 (2013).</footnote>

<footnote n="67">⁶⁷M. Ruckenbauer, S. Mai, P. Marquetand, and L. González, "Revealing deactivation pathways hidden in time-resolved photoelectron spectra," *Sci. Rep.* **6**, 35522 (2016).</footnote>

<footnote n="68">⁶⁸S. Mai, M. Richter, P. Marquetand, and L. González, "The DNA nucleobase thymine in motion—Intersystem crossing simulated with surface hopping," *Chem. Phys.* **482**, 9–15 (2017).</footnote>

<footnote n="69">⁶⁹T. J. Martínez, "Ab initio molecular dynamics around a conical intersection: Li(2p) + H₂," *Chem. Phys. Lett.* **272**, 139–147 (1997).</footnote>

<footnote n="70">⁷⁰I. Fdez. Galván, M. G. Delcey, T. B. Pedersen, F. Aquilante, and R. Lindh, "Analytical state-average complete-active-space self-consistent field nonadiabatic coupling vectors: Implementation with density-fitted two-electron integrals and application to conical intersections," *J. Chem. Theory Comput.* **12**, 3636–3653 (2016).</footnote>

<footnote n="71">⁷¹M. A. El-Sayed, "Spin-orbit coupling and the radiationless processes in nitrogen heterocyclics," *J. Chem. Phys.* **38**, 2834–2838 (1963).</footnote>

<footnote n="72">⁷²M. Baba, "Intersystem crossing in the $1n\pi^*$ and $1\pi\pi^*$ states," *J. Phys. Chem. A* **115**, 9514–9519 (2011).</footnote>
<footnote n="48">⁴⁸F. Aquilante, T. B. Pedersen, and R. Lindh, "Low-cost evaluation of the exchange Fock matrix from Cholesky and density fitting representations of the electron repulsion integrals," *J. Chem. Phys.* **126**, 194106 (2007).</footnote>

<footnote n="49">⁴⁹J. Boström, M. G. Delcey, F. Aquilante, L. Serrano-Andrés, T. B. Pedersen, and R. Lindh, "Calibration of Cholesky auxiliary basis sets for multiconfigurational perturbation theory calculations of excitation energies," *J. Chem. Theory Comput.* **6**, 747–754 (2010).</footnote>

<footnote n="50">⁵⁰B. A. Hess, C. M. Marian, U. Wahlgren, and O. Gropen, "A mean-field spin-orbit method applicable to correlated wavefunctions," *Chem. Phys. Lett.* **251**, 365–371 (1996).</footnote>

<footnote n="51">⁵¹P.-À. Malmqvist and B. O. Roos, "The CASSCF state interaction method," *Chem. Phys. Lett.* **155**, 189–194 (1989).</footnote>

<footnote n="52">⁵²P.-À. Malmqvist, B. O. Roos, and B. Schimmelpfennig, "The restricted active space (RAS) state interaction approach with spin-orbit coupling," *Chem. Phys. Lett.* **357**, 230–240 (2002).</footnote>

<footnote n="53">⁵³G. Kovacevic and V. Veryazov, "Luscus: Molecular viewer and editor for MOLCAS," *J. Cheminf.* **7**, 16 (2015).</footnote>

<footnote n="54">⁵⁴T. J. Martínez, M. Ben-Nun, and R. D. Levine, "Molecular collision dynamics on several electronic states," *J. Phys. Chem. A* **101**, 6389–6402 (1997).</footnote>

<footnote n="55">⁵⁵B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martínez, "Implementation of ab initio multiple spawning in the Molpro quantum chemistry package," *Chem. Phys.* **347**, 3–16 (2008).</footnote>

<footnote n="56">⁵⁶J. D. Coe, B. G. Levine, and T. J. Martínez, "Ab initio molecular dynamics of excited-state intramolecular proton transfer using multireference perturbation theory," *J. Phys. Chem. A* **111**, 11302–11310 (2007).</footnote>

<footnote n="57">⁵⁷H. Tao, B. G. Levine, and T. J. Martínez, "Ab initio multiple spawning dynamics using multi-state second-order perturbation theory," *J. Phys. Chem. A* **113**, 13656–13662 (2009).</footnote>

<footnote n="58">⁵⁸H. Tao, T. K. Allison, T. W. Wright, A. M. Stooke, C. Khurmi, J. van Tilborg, Y. Liu, R. W. Falcone, A. Belkacem, and T. J. Martinez, "Ultrafast internal conversion in ethylene. I. The excited state lifetime," *J. Chem. Phys.* **134**, 244306 (2011).</footnote>

<footnote n="59">⁵⁹T. Mori, W. J. Glover, M. S. Schuurman, and T. J. Martinez, "Role of Rydberg states in the photochemical dynamics of ethylene," *J. Phys. Chem. A* **116**, 2808–2818 (2012).</footnote>

<footnote n="60">⁶⁰T. S. Kuhlman, W. J. Glover, T. Mori, K. B. Møller, and T. J. Martínez, "Between ethylene and polyenes—the non-adiabatic dynamics of cis-dienes," *Faraday Discuss.* **157**, 193–212 (2012).</footnote>

<footnote n="61">⁶¹T. J. A. Wolf, T. S. Kuhlman, O. Schalk, T. J. Martínez, K. B. Møller, A. Stolow, and A.-N. Unterreiner, "Hexamethylcyclopentadiene: Time-resolved photoelectron spectroscopy and ab initio multiple spawning simulations," *Phys. Chem. Chem. Phys.* **16**, 11770–11779 (2014).</footnote>

<footnote n="62">⁶²L. Liu, J. Liu, and T. J. Martinez, "Dynamical correlation effects on photoisomerization: Ab initio multiple spawning dynamics with MS-CASPT2 for a model trans-protonated Schiff base," *J. Phys. Chem. B* **120**, 1940–1949 (2016).</footnote>

<footnote n="63">⁶³W. J. Glover, T. Mori, M. S. Schuurman, A. E. Boguslavskiy, O. Schalk, A. Stolow, and T. J. Martínez, "Excited state non-adiabatic dynamics of the smallest polyene, trans 1,3-butadiene. II. Ab initio multiple spawning simulations," *J. Chem. Phys.* **148**, 164303 (2018).</footnote>

<footnote n="64">⁶⁴J. W. Park and T. Shiozaki, "On-the-fly CASPT2 surface-hopping dynamics," *J. Chem. Theory Comput.* **13**, 3676–3683 (2017).</footnote>

<footnote n="65">⁶⁵M. Richter, S. Mai, P. Marquetand, and L. González, "Ultrafast intersystem crossing dynamics in uracil unravelled by ab initio molecular dynamics," *Phys. Chem. Chem. Phys.* **16**, 24423–24436 (2014).</footnote>

<footnote n="66">⁶⁶S. Mai, P. Marquetand, M. Richter, J. González-Vázquez, and L. González, "Singlet and triplet excited-state dynamics study of the keto and enol tautomers of cytosine," *ChemPhysChem* **14**, 2920–2931 (2013).</footnote>

<footnote n="67">⁶⁷M. Ruckenbauer, S. Mai, P. Marquetand, and L. González, "Revealing deactivation pathways hidden in time-resolved photoelectron spectra," *Sci. Rep.* **6**, 35522 (2016).</footnote>

<footnote n="68">⁶⁸S. Mai, M. Richter, P. Marquetand, and L. González, "The DNA nucleobase thymine in motion—Intersystem crossing simulated with surface hopping," *Chem. Phys.* **482**, 9–15 (2017).</footnote>

<footnote n="69">⁶⁹T. J. Martínez, "Ab initio molecular dynamics around a conical intersection: Li(2p) + H₂," *Chem. Phys. Lett.* **272**, 139–147 (1997).</footnote>

<footnote n="70">⁷⁰I. Fdez. Galván, M. G. Delcey, T. B. Pedersen, F. Aquilante, and R. Lindh, "Analytical state-average complete-active-space self-consistent field nonadiabatic coupling vectors: Implementation with density-fitted two-electron integrals and application to conical intersections," *J. Chem. Theory Comput.* **12**, 3636–3653 (2016).</footnote>

<footnote n="71">⁷¹M. A. El-Sayed, "Spin-orbit coupling and the radiationless processes in nitrogen heterocyclics," *J. Chem. Phys.* **38**, 2834–2838 (1963).</footnote>

<footnote n="72">⁷²M. Baba, "Intersystem crossing in the $1n\pi^*$ and $1\pi\pi^*$ states," *J. Phys. Chem. A* **115**, 9514–9519 (2011).</footnote>