Electronic Supplementary Information Nonadiabatic dynamics in multidimensional complex potential energy surfaces

S3 Computational details for the application to iodoethene 3 S3.1 Bound state calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 S3.2 Scattering calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 S3.3 Dynamics simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 S3.4 DEA cross section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 S3.5 Modelling the complex PES of iodoethene . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5


S1 Generalised surface hopping method
For an arbitrary set {ψ j } which is used to expand the electronic wave function (eq. 2 in the main text), the fewest switches probabilities 1 are generalised to P j→k = max 2∆t When the set {ψ j } corresponds to the eigenstates of the real-valued electronic Hamiltonian H R , then the dynamics takes place in the adiabatic representation. In this case, the adiabatic potential energy surfaces (PESs) are coupled by the nonadiabatic coupling vectors F jk . Alternatively, one could work in a basis where the nonadiabatic couplings are minimised, but then the diabatic PESs would interact via the off-diagonal diabatic couplings H jk . The adiabatic representation tends to be preferred in on-the-fly trajectory methods, since the nonadiabatic couplings are more strongly peaked when compared to the smoother diabatic couplings 2 . On top of that, the adiabatic energies and nuclear gradients are directly provided when solving the electronic problem, dispensing diabatisation techniques. Therefore, in the adiabatic representation, H R jk = 0 for j = k, and eq. S1 reduces to eq. 5 of the main text.
Concerning the decoherence correction, we have employed a modified version of the simplified decay of mixing 3 . At each time step, the coefficients c were updated according to where j is the current state, E kin is the nuclear kinetic energy, and α is a parameter that controls the decoherence time τ k j . The difference to the usual simplified decay of mixing 3 lies in eq. S3. After damping the other coefficients (eq. S2), c j is renormalised according to the total population p = ∑ k |c k | 2 , which in general is different from unity. Once the dynamics simulation is finished, any time-dependent observable of interestŝ(t) can be computed by performing a weighted average over the contribution from each trajectory, aŝ where the sum runs over a total of N t trajectories, s i (t) represents the observable and p i (t) the total population, for trajectory i and at time t.

S2 Propagation on the model potential energy surfaces
The following set of parameters was adopted for the 2REH model (defined in eqs. 6-11 of the main text): , and a reduced mass of µ = 9 amu. Weaker (V 12 (R c ) = 0.03 eV) and stronger (V 12 (R c ) = 0.1 eV) diabatic coupling strengths have been considered.
The adiabatic representation was employed for the propagation with the CS-FSSH method. For each initial resonant state, 1,000 trajectories were propagated, with a time step of 0.1 fs for the classical part, and of 0.005 fs for the TDSE. Initial conditions were generated by sampling from the Wigner distribution of the vibrational ground state of V 0 . The fifth order Butcher's integrator 4 was employed for integrating of the TDSE, with a time step of 0.005 fs. The classical integration was performed with the velocity Verlet algorithm, with a time step of 0.1 fs. For the shorter steps, energies, nuclear velocities, nonadiabatic couplings, and resonance widths were obtained by interpolating the values computed at the longer steps of the classical propagation. After a successful hopping, the momentum was adjusted in order to conserve total energy, and for a frustrated hopping, the momentum was maintained. A decoherence correction with α = 0.1 Hartree 5 was employed.
The quantum wavepacket was propagated on the diabatic PESs, as described elsewhere 6,7 . The propagation was performed with the split-operator technique 8 , with a spacial grid of 2 14 points, from −2 a 0 to 14 a 0 , and a time grid of 2 18 points, from 0 to 500 fs.

S3 Computational details for the application to iodoethene
The description of transient anion states should be performed with either scattering methods [9][10][11] or adapted bound state methods [12][13][14] . Alternatively, usual bound state (quantum chemistry) calculations can also provide a quantitative description of such states [15][16][17][18][19][20] . The latter option must be used with caution, though, as pseudo-continuum states associated with the free electron could appear as a solution of the electronic Hamiltonian. In the present application, we have employed a mixed approach. Conventional bound state methods have been employed to compute resonance energies, albeit with a correction from scattering calculations.

S3.1 Bound state calculations
During the dynamics propagation, anion and neutral energies were computed with the multireference configuration interaction (MRCI) level of theory, with the latest implementation of the Columbus package 21 . The active space comprised 6 orbitals (the bonding σ CI , two bonding π, the non-bonding n I , and the antibonding π * and σ * CI orbitals), which shared 8 electrons in the neutral and 9 in the anion. The reference wavefunction was obtained with a state-average complete active space self-consistent field calculation (SA-CASSCF) for the three lower-lying neutral singlet states in the state averaging, while considering the same (8,6) active space. We have employed the cc-pVDZ basis set 22 for carbon and hydrogen atoms, and the jorge-ADZP basis set 23 for the iodine atom, as implemented in Columbus 21,24 . The MRCI wavefunction of both anion and neutral states were obtained by accounting for single excitations from the doubly occupied orbitals (excluding the core) into the active space, and from the latter into the all virtual orbitals.

S3.2 Scattering calculations
The fixed-nuclei elastic scattering calculations were performed with the Schwinger multichannel method (SMC) 11 . The scattering wavefunction was expanded in a set of configuration state functions (CSFs), which were represented as an antisymmetrised product of a target function and a scattering orbital. Here, we have considered CSFs built from the neutral ground state (as described in the restricted Hartree-Fock approximation), and from singlet and triplet singly excited functions. Modified virtual orbitals 25 were employed to represent the unoccupied orbitals, as obtained in the field of a +8 cationic Fock operator. The CSFs were selected according to an energy criterion for the relevant orbitals 26 , with a cutoff energy of 1.5 Hartree. Pseudopotentials for the iodine and carbon atoms were employed for the core electrons 27 , while valence electrons were described with the same Cartesian Gaussian functions as presented elsewhere 28 . This set has been supplemented with 3s and 3p functions centred at carbon bonded to the iodine atom, with eventempered exponents generated with a ratio of 4.

S3.3 Dynamics simulations
A development version of the Newton-X package 29,30 was employed for the dynamics simulations. In total, 1,600 trajectories were propagated, each half starting at one anion state. The initial conditions were generated by sampling from the Wigner distribution of the quantum harmonic oscillator associated with each vibrational normal mode 2,31 , for a temperature of 333 K. Reasonably converged observables were attained when sampling from the Wigner distribution, differently from a previous dynamics simulation of DEA 20 , where modified sampling techniques 32 were required. Newton's equations of motion were integrated with the velocity Verlet algorithm, with a time step of 0.25 fs. At each step, an additional energy calculation for the neutral ground state was performed as well, which was required to evaluate the resonance widths, as detailed in Sec. S3.5. Meanwhile, the electronic TDSE was propagated according to the fifth-order Butcher's integrator 4 , with a step of 0.002 fs. For the shorter steps, energies, nuclear velocities, nonadiabatic couplings, and resonance widths were obtained by interpolating the values computed at the longer steps of the classical propagation. The time step for the TDSE is one order of magnitude shorter than what is usually employed for dynamics on real PESs, which was found to be critical for properly describing the effect of the imaginary component in the present case. When hopping between surfaces took place, nuclear momenta were adjusted in the direction of the nonadiabatic coupling vector 33 , and in the few trajectories that incurred into frustrated hopping, they were kept the same. The standard value of α = 0.1 Hartree 5 was employed in the expression for the decoherence time (eq. S4). The propagation was terminated when the C−I distance reached 6a 0 , as long as at least 25 fs had passed.

S3.4 DEA cross section
In our application, we are particularly interested in computing the DEA cross section σ k associated to each precursor anion state k, which can be computed 20,34 (in atomic units) according to In the above expression, E is the electron impact energy, p k,i (t → ∞) is the final population of trajectory i, and the first term in the summation denotes an energy-dependent resonance width, which is proportional to the probability of electron capture. It is comprised of the resonance width at t = 0, Γ 0 k,i , and an energydependent component γ k (E). The latter depends on the collision energy E, in the numerator, and on the resonance energy at t = 0, E 0 k,i , in the denominator. These functions were obtained by fitting the equilibriumgeometry eigenphase sum (obtained with the scattering calculations), according to 35 with where erf is the error function. The set of obtained parameters for the lower (k = 1) and upper (k = 1) anion states are shown in Tab. S1. Finally, g k,i (E) was represented by a normalised Gaussian function centred at the resonance energy E 0 k,i , with a line width of 0.05 eV.  Table S1 Parameters of the energy-dependent functions γ 1 (E) and γ 2 (E) (for E given in eV), which appear in the expression for the DEA cross sections (eq. S6). Fig. S1 shows the set of results that were required in order to build the autodetachment model for the dynamics. First, it shows the elastic cross sections obtained at the neutral equilibrium geometry, as computed with the SMC scattering calculations. The peaks centred at 0.58 eV and 1.11 eV correspond to the σ * and π * shape resonances, respectively, which match quite well with electron transmission spectra results (0.5 eV and 1.05 eV) 36 . The scattering calculations also provide vertical resonance widths of 0.16 eV for the σ * state and of 0.26 eV for the π * state, corresponding to autodetachment lifetimes of 4.2 fs and 2.5 fs.  Figure S1 On the top left panel, symmetry resolved elastic cross sections, computed at the neutral equilibrium geometry. On the top right, behaviour of the model resonance widths for the lower Γ 1 and higher-lying Γ 2 anion states, as described in the text. On the bottom, PESs along the C−I and C=C stretching from the neutral equilibrium values, for the neutral, σ * and π * energies (triangles) and resonance widths obtained from the scattering calculations (circles) and from the model Γ 1 and Γ 2 (dashed curves).

S3.5 Modelling the complex PES of iodoethene
Despite the good agreement between SMC calculations and experiment, the resonance energies computed with the MRCI method are expected to be overestimated, given the well-known 10,14 difficulty in obtaining an accurate and balanced description of correlation effects for anion and neutral states. Indeed, we have obtained 2.11 eV and 3.16 eV for the σ * and π * states at the neutral equilibrium geometry, compared to the experimentally observed 0.5 eV and 1.05 eV 36 . Due to inaccuracies on the employed electronic structure level of theory, systematic shifts or scaling of the computed PESs are sometimes performed 20,37,38 , and a similar procedure was adopted here. The resonance energy E j associated to the anion state j was evaluated is the MRCI resonance energy for the anion state j, given by the difference between anion V j and neutral V 0 energies, as computed along the dynamics, and ∆E is a constant energy shift. This correction is such that the vertical (R = 0) resonance energy E j for the state that triggers the dynamics j matches the value obtained with the ab-initio SMC calculation E SMC j (0). Accordingly, when the dynamics started from the lower-lying anion state, the correction was computed as (0), and when the higher-lying state was initially formed, as ∆E = E SMC 2 (0) − E MRCI 2 (0). We are aware that such an unbalance between neutral and anion states at the MRCI level is excessive, but as our main goal in this paper is to demonstrate the potential of the CS-FSSH methodology for dynamics simulations of real molecular systems. Moreover, we recall that this unbalance comprises a well-recognised and present challenge for electronic structure theory 10,14 . One has to move to very accurate descriptions of electronic correlation in order to achieve reasonably converged resonance energies, by accounting for triple excitations in coupled cluster methods for instance 39 . In the present case, when we include double excitations as well as the Davidson correction 40 (MRCISD+Q), the average deviation to experiment drops from 1.9 eV to 0.7 eV. Discrepancies around 0.5 eV have already been observed in calculations of resonant anions employing the MRCISD+Q 20 and the CASPT2 methodology [15][16][17] . Unfortunately, the MRCISD+Q level of theory does not count with analytical nuclear gradients, which is a required ingredient for the nuclei propagation. Furthermore, as our main goal in this paper is to demonstrate the use of the CS-FSSH methodology, we are satisfied that the PES shift procedure can provide a fair description of the resonances of iodoethene.
For the dynamics simulations, the imaginary PESs were modelled based on the following series of arguments. Since resonance energies and resonance widths are usually positively correlated, for both σ * 20,41,42 and π * 43 states, here we have assumed that the resonance widths of each state depend only on the corresponding resonance energy, i.e., Γ j (R) Γ j [E j (R)]. Now, the resonance widths should be described more accurately in regions of the PESs where they are larger, because there the population decays more rapidly. Once the resonant anion state is formed, its energy tends to decrease, hence also the resonance energy and width. Therefore, it is more advisable to have an accurate model for the resonance widths at the beginning of the dynamics, when they are larger; than at later stages, when they would be smaller and thus of diminishing relevance to the population decay.
From a small sample of very short trajectories, the σ * resonance was found to stabilise by C−I bond stretching, while the π * resonance relaxed by stretching of the C=C bond. That being so, we have performed a series of scattering calculations along both stretching coordinates (the anion PESs are shown in Fig. S1), which served to build an approximated model for the widths of the lower Γ 1 and upper Γ 2 anion states, as follows. The width of the lower-lying state, Γ 1 , was obtained from the behaviour of the σ * widths along the C−I coordinate (dashed blue curve in the bottom left panel of Fig. S1). Similarly, the behaviour of the π * width along the C=C coordinate (dashed red curve in the bottom right panel) provided the width of the higher-lying state, Γ 2 . More precisely, the instantaneous resonance widths adopted for the dynamics were modelled with a third-order polynomial interpolation 44 to these pre-computed pairs of energies and widths, The set of polynomial parameter employed for each energy interval is presented in Tab. S2, and the corresponding Γ 1 (E 1 ) and Γ 2 (E 2 ) are shown in the top right panel of Fig. S1.  Table S2 Parameters of the polynomial interpolation (eq. S10), for energies and widths given in eV.
As additional validation of our assumptions, the profile of each Γ was checked against the coordinate that was not employed to build it in the first place. In the bottom panels of Fig. S1, dashed curves represent the model Γ 1 and Γ 2 , obtained as described above, while dots are the data from the scattering calculations. We have found a good matching between Γ 1 and the σ * widths along the C=C stretching, and also between Γ 2 and the π * widths along the C−I stretching. This result further supports the above assumptions for the resonance widths and suggests that the conceived autodetachment model should be accurate enough for our purposes.