Probing deactivation pathways of DNA nucleobases by two-dimensional electronic spectroscopy: first principles simulations.

The SOS//QM/MM [Rivalta et al., Int. J. Quant. Chem. 2014, 114 , 85.] method consists of an arsenal of computational tools allowing accurate simulation of onedimensional (1D) and bidimensional (2D) electronic spectra of monomeric and dimeric systems with unprecedented details and accuracy. Prominent features like doubly excited local and excimer states, accessible in multi-photon processes, as well as charge-transfer excimer states arise naturally through the fully quantum-mechanical description of the aggregates. In this contribution the SOS//QM/MM approach is extended to simulate time-resolved 2D spectra that can be used to characterize ultrafast excited state relaxation dynamics with atomistic details. We demonstrate how critical structures on the excited state potential energy surface, obtained through state-of-the-art quantum chemical computations, can be used as snapshots of the excited state relaxation dynamics to generate spectral ﬁngerprints for different de-excitation channels. The approach is based on high-level multi-conﬁgurational wavefunction methods combined with non-linear response theory and incorporates the effects of the sol-vent/environment through hybrid quantum mechanics/molecular mechanics (QM/MM) techniques. Speciﬁcally, the protocol makes use of the second-order Perturbation theory (CASPT2) on top of Complete Active Space Self Consistent Field (CASSCF) strategy to compute the high-lying


Introduction
Ultrafast pump-probe spectroscopy experiments backed up by theoretical models at quantum-mechanical level have achieved a detailed understanding of the photophysics of the DNA building blocks adenine, guanine, thymine and cytosine in gas-phase and solution. A thorough overview of the current status of research is given in several review articles [1][2][3][4] . Multimeric nucleobase systems formed through stacking (single-stranded DNA) and base-pairing (double-standed DNA) add many new facades to the photophysics of single bases: formation and dynamics of exciton states 5,6 , population and decay of intra-and inter-strand charge transfer (CT) excimer states 7,8 , charge separation and recombination [9][10][11] , ultrafast decay through photoinduced proton transfer 7,[12][13][14][15] , formation of lesions [16][17][18] . Interpretation of the electronic spectra of DNA multimers, which is the entirety of all these processes, is a challenge for the conventional one-dimensional (1D) pump-probe spectroscopy. When multiple electronic transitions come into resonance with the excitation wavelength the spectroscopic signatures of individual de-excitation pathways overlap making their disentanglement virtually impossible. Population transfer between different channels cannot be tracked because the correlation between pump and probe wavelengths is not resolved. This complexity calls for a spectroscopic method with an enhanced spectral resolution. Ultrafast bidimensional (2D) spectroscopy holds the premise to disentangle multiple deactivation pathways. First applications in the visible (Vis) and in the near-UV (NUV) have been documented, resolving energy transfer mechanisms in light harvesting proteins 19,20 and solvation dynamics in nucleotides 21,22 . The evolution of experimental techniques raises the expectations to theoretical modeling of nonlinear spectra. The theoretical framework of nonlinear electronic spectroscopy has been worked out [23][24][25] and nowadays methods are sought to feed the resulting equations with system-specific data [26][27][28] . Static and dynamic studies at the quantum-mechanics level have been successfully applied to study the excited state potential energy surfaces of DNA dimers in gas-phase 8 , in aqueous solution 29 and in the native environment of the solvated single-and doublestranded DNA [30][31][32] , thus providing a tool for obtaining information about excited state energetics and geometrical deformations. Time-resolved Vis-pump IR/Visprobe 33 and UV-pump stimulated-Raman-probe 28 spectra simulations demonstrate how quantum-mechanical data can be translated into experimental observables. Simulating ultrafast 2D electronic spectroscopy from first principles is challenging as it requires to incorporate knowledge about the relaxation dynamics of singly and doubly excited states accessible along the de-excitation pathways into the non-linear spectroscopy simulations. Then, it is the excited state absorption (ESA) and the stimulated emission (SE) whose temporal evolution re-flects the geometrical and electronic changes along the potential energy surface. Novel computational and conceptual strategies are required. The SOS//QM/MM protocol is the first method which couples quantum-mechanics / molecular mechanics (QM/MM) calculations, multi-configurational excited state calculations with non-linear spectroscopy to generate time-resolved 2D electronic spectra 34 . It allows for the characterization of the vertical 1D and 2D electronic spectra of monomeric and dimeric systems with unprecedented details and accuracy 27,35,36 . Prominent features like doubly excited local and excimer states, accessible in multi-photon processes, as well as CT excimer states are revealed. In this paper we present an extension of this protocol to simulate coherent excited state dynamics based on static and dynamic computational data 37 . Under the assumption that the system's time evolution is slow compared to signal generation, we are able to probe single snapshots along the excited state relaxation pathways, which can be achieved experimentally in a pump-probe pulse arrangement through scanning the delay time between the pump and the probe pulse sequences. The protocol is used to generate the spectroscopic signatures at selected geometries along the major deactivation channels of the stacked adenine-adenine (A-A) dimer in an explicitly solvated double-stranded DNA model.  We treat the molecular system as a three-level model, whose eigenstates can be formally devided into: g the ground state; e a set of singly excited states, accessible via the pump pulse sequence; f a set of higher excited states, which are probed after a time delay t 2 . The nonlinear response functions for ground state bleaching (GSB), SE and ESA, whose Feynman pathways are shown in Fig. 1 for the rephasing pulse orientation technique, read 23,24,38 : with Θ(t) the Heavyside step-function, µ ab the vector of the transition dipole moment between states a and b, ω ab the frequency of the electronic transition between states a and b, γ a phenomenological dephasing rate constant and G e→e (t 2 ) a function describing the evolution of the system during the time interval t 2 . The transition e → e could describe both evolution in the same electronic state and population transfer to another electronic state. It determines the electronic spectrum of the system after the delay time t 2 which enters eq. 1 through transition frequencies and dipole moments (shown in red). Ideally, quantum or mixed quantum-classical dynamics simulations are conducted to obtain these parameters. Parametrized models which solve the secular Redfield equations have been used to account for population relaxation and Stokes shifts at early times t 2 without explicit dynamics simulations 39 . The temporal evolution of the nonlinear response is a function of the electronic and geometrical changes in the system modulated by the surrounding. Due to their nature GSB signals do not exhibit time-dependent spectral shifts and their recovery is determined by the decay rates of all excitations resonant with the pump pulse. These features make bleach signals excellent probes for extracting decay rates by means of time-resolved 1D ultrafast spectroscopy. However, the correlation between decay rates and de-excitation pathways cannot be resolved. On the other hand, SE and ESA are specific to each excited state, which makes them eligible probes for pathway-specific extraction of decay rates. Tracking SE and ESA signals by means of 1D ultrafast spectroscopy is challenging as signals from competing channels overlap spectrally and temporally. Furthermore, the temporal evolution of the SE and ESA often span over a large spectral range. 2D time-resolved spectroscopy utilizing broadband probe pulses disentangles spectrally and resolves temporally the fingerprints of the individual de-excitation channels [40][41][42] . In this work we employ geometry optimizations in the excited state and probe the higher excited manifold f and the ground state g at the obtained stationary points, while keeping the pump pulse pair in resonance with the singly excited state manifold e at the Franck-Condon (FC) point. The local minima are believed to trap excited state population and are, thus, responsible for decay times in the picosecond (ps) time range. With this approach the characteristic spectral signatures of different de-activation channels can be obtained. To compute the fingerprints of each de-excitation channel we regard the channels as independent, i.e. we assume that no population is transferred between them or to the ground state (GS) during relaxation from the FC region. With this simplification the 2D spectrum of the FC point can be formally divided along ω 1 into traces resonant with either transition. The time-dependent variation of the traces is analyzed by comparing the signals at the FC point and at the stationary points obtained in this study. The following approximations are adopted: the system evolves coherently along the populated excited state, i.e. inhomogeneous broadening is neglected; coherence dynamics during t 1 and t 3 is neglected, i.e. the system's time evolutuon is assumed to be slow compared to the signal generation; inter-state coherences (i.e. the system is not in a population state during t 2 in the SE and ESA diagrams, Fig. 1) are assumed to decay prior to probing. Note that geometry optimizations do not contain temporal information and, therefore, we cannot provide time delays. For simulating quasi-absorptive (rephasing and non-rephasing) 2D electronic spectra we have combined the QM/MM methodology, which provides the transition dipole moments at the Complete Active Space Self Consistent Field (CASSCF) level 43 and transition energies corrected at the second-order Perturbation theory (CASPT2) 44 , with the sum-over-states (SOS) approach 45 . SOS calculations were performed with Spectron 2.7 24 readapting the energy levels calculated at each stationary point in order to include the GS bleaching at the FC, i.e. by maintaining the FC g-e energy gap and rigidly shifting all the f energies calculated at the stationary point in order to align the e energy with the FC e energy. Using the fingerprint traces and adopting some assumptions supported by experimental data and the mechanistic picture obtained at the ab initio computational level ideal 2D spectra for different scenarios can be simulated. These spectra show how static data from excited state optimizations and minimum energy path calculations can be translated into spectroscopic signatures.
3 Application to the adenine-adenine dimer in solvated DNA 3.1 Experimental and theoretical observations for poly(dA) systems Femtosecond (fs) transient absorption experiments by Kohler and co-workers support a model in which DNA multimer excitations decay to long lived excimer states with a decay time constant ranging from ten to hundreds of ps 46 . Adenine homopolymers and homoduplexes made of adenine-thymine base pairs show much slower decay times compared to those of isolated bases 1,2,47,48 . More recent findings support the fact that long lived excited states decay more rapidly in double-stranded poly(dA-dT) than in single stranded poly(dA) sequences 49,50 , addressing the role of inter-strand proton transfer in the deactivation mechanism 7,12,13 . Long-lived signals have also been probed by means of fluorescence up-conversion techniques 5,51,52 , and more recently through femtosecond infrared spectroscopy by probing the distinct cationic and anionic species formed in the process 10 . The experimental findings suggest that excited states localized on just two stacked bases are the common trap states independently of the number of stacked nucleotides 6,53 . On the other hand, Markovitsi and co-workers have postulated the possible delocalization of the exciton over several nucleobases 5,51,52 (up to six according to some theoretical models 54 ) by means of fluorescence upconversion experiments. They suggest that the long-lived component is caused by charge relocalization on the excited monomers. The de-excitation channels present in systems such as the adenine-thymine base pair have been recently studied with 2D photon echo IR spectroscopy, yielding novel fingerprints that embody the different motions undergone by the molecules during the relaxation process 55,56 . Several computational works have been undertaken on DNA multimer models to elucidate how light interacts with these electronically and structurally complex systems. Merchán and co-workers computed energies and structure of adenine homodimers in vacuum and in water at CASPT2 level 8 . They support the idea that the long-lived excited states seen in DNA are of excimer type, formed by stacking of two bases. Two stacked conformations wre characterized, a perfectly stacked "face-to-face" orientation and a poorly stacked orientation where both nucleobases are twisted. Going from the poorly stacked towards the perfectly stacked conformation the CT states decrease in energy, therefore concluding that the population of the CT states upon photoexcitation strongly depends on the conformational properties of DNA. They associate the sub-ps life-time with two parallel routes: one leads to the formation of the excimer states, the other to the decay to the ground state directly from L a state. The formed excimer state would display much slower decay times (i.e. in the ps regime), as the system has to overcome an energy barrier to access a conical intersection (CI) with the GS. A time-dependent DFT (TD-DFT) study by Barone and co-workers concerning different single-and double-stranded poly(dA) multimers embedded in water (treated as a polarizable continuum) also indicate that the long living component of the excited state population correspond to a dark excimer produced by a charge transfer between stacked adenines 29 . In their study the CT excimer constitutes the absolute excited state minimum. A QM/MM study by Plasser et al. utilizing the multi-reference configuration interaction (MRCI) ab initio method on the adenine dinucleotide in water predicts the formation of several stable low lying exciplexes of nπ * and ππ * character with short inter-molecular separation but without notable charge transfer. The authors suggest that the long life-time results from trapping the system in these minima 31 . Based on QM/MM computations of the adenine monomer and homodimer embedded in a double-stranded DNA, instead, Conti et al. propose that an intra-monomer mechanism associated with relaxation of the L a channel is compatible with the observed multi-exponential decay, including the longer (> 100 ps) life-time component 57,58 . Summing up, there is still no agreement on the de-excitation mechanisms in nucleobase multimers, particularly regarding the participation of CT excimer states.

Mechanistic picture and de-excitation pathways
In the following we outline the recent results of Conti et al. for the adenine dimer embedded in a double-stranded DNA model 58 . Both nucleobases were treated quantum-mechanically at CASPT2//CASSCF level with the surrounding DNA chain and solvent represented through point charges. The bright excited states with L a (3') and L a (5') character (3' and 5' denote the upper and lower nucleobase in the stacked dimer in Fig. 2) couple and form an exction pair L a (3')-L a (5')/L a (3')+L a (5') which absorbs around 5.00 eV with a small splitting of 0.05 eV (500 cm −1 ). After excitation the exciton states evolve toward monomer localized excitations with a deactivation mechanism similar to the one reported in gas-phase, occurring through a C 2 puckering distortion 59-63 (see Fig. 2 ∼65°at the CI). After the initial energy minimization, associated with bond rearrangements, a broad plateau is reached (Fig. 2). Only at large twisting angles (i.e. above 50°) the excited state energy decreases again towards the CI. Conti et al. suggest that the lack of a steep gradient towards the CI may be the cause for the long life-time encountered in poly(dA) multimers.
At the relaxed GS geometry the lowest pair of CT states (H(3')→L(5') and deformation. Both CIs are accessible only via barriers which offer another potential source of the long life-time encountered in poly(dA) multimers 5,46,51,52,65 . We note, however, that the barriers evaluated by Conti et al. are too large (∼ 1 eV) to correlate with ps life-times. Furthermore, despite a thorough search conformations with a CT state below 5.50 eV could not be located.
In the course of geometry optimization Conti et al. observe that the C 2 puckering distortion is a common feature of both the L a and CT relaxation pathways, even though the C 2 puckering along the CT pathway is less pronounced. Along the L a deactivation pathway the puckering is a result of the ethylene-like twist arond the C 2 -C 3 double bond, necessary to reach the CI. Along the CT deactivation pathway the puckering reduces the electron repulsion of the excess charge density in the negatively charged monomer. A 2D map in the space of the puckering and inter-base distance coordinates indicates that population transfer through vibrational energy redistribution can occur in both directions via small barriers 58 . The oscillation between deactivation channels is likely to slow down the de-excitation process.

Computational details
The critical geometries obtained in the mechanistic study of the bright exciton states L a (3')-L a (5') and L a (3')+L a (5') and of the lowest CT states CT ( Calibration against full active space calculations of the adenine monomer in gasphase shows that this highly reduced AS does not permit quantitative predictions (details in the Supplementary Information). In fact, a number of transitions accessible by visible light out of the L a band, which involve lower occupied and higher virtual π-orbitals, are missed. Despite this, the reference calculations demonstrate that doubly excited states which arise via excitations among H-1, H, L and L+1 acquire the highest oscillator strengths and, thus, dominate the NUV spectral window. We expect these excited states to be affected the most along the de-excitation pathways, considering that they involve the same orbitals as the optimized L a and CT states. Thus, our efforts focus in resolving the qualitative spectral dynamics of the ESA into the doubly excited states in the NUV spectral window between 20000 cm −1 and 40000 cm −1 from the L a band (i.e. between 60000 cm −1 and 80000 cm −1 from the GS). A shift of 3000 cm −1 was applied to all states (implies a 6000 cm −1 shift of all mixed doubly excited states) in order to better reproduce the reference calculations. This shift has no effect on the position of the SE and ESA in the spectra, but blue-shifts the GSB by 3000 cm −1 .
The generally contracted ANO-L basis set was utilized and the following contraction scheme was adopted:

Traces for L a and CT
First, we discuss the traces of the individual excited state relaxation pathways out of the FC region towards the local minima.
We begin with the trace of the L a de-excitation channel from the FC region to To deal with this issue we selected two representative geometries with C 2 puckering 11°and 40°for generating the spectra. The L a band gives rise to the most intense transitions in the NUV, absorbing around 39000 cm −1 . They are delocalized over both nucleobases due to exciton coupling forming the exciton states L a (3')-L a (5') and L a (3')+L a (5'). Fig. 3 shows the positive linear combination which collects most of the oscillator strength in the NUV. The region around the GSB of the L a (3')+L a (5') (37000-41000 cm −1 ) is congested with off-diagonal bleach signals (e.g. peak 2) and ESA to mixed doubly excited states (e.g. L a (3'+5') peak 5) and is hard to disentangle. Therefore, we encourage to perform experiments in the low energy window below the bleach region. Pairs of ESA peaks appear (peaks 3|4, and 6|7), which 1-18 | 9 arise via probing the local doubly excited states D * 1 (3'/5') and D * 2 (3'/5') out of the L a band. The ESA doublets 3|4, and 6|7 split as the system leaves the FC region and the exciton localizes on one of the bases. Thereby, the ESA to doubly excited states on the puckered adenine exhibit no shift (peak 3) or red-shifts (peaks 6 and 8) while ESA to doubly excited states on the unmodified adenine exhibit a blue-shift (peaks 4 and 7) equal to the red-shift of the SE (peak 1'). The pronounced red-shift of peaks 6 and 8 indicates that the geometrical deformations along the L a relaxation pathway stabilize local doubly excited states even stronger which is a remarkable finding. The SE red-shifts continuously even as the excited state reaches a plateau since the GS exhibits constant destabilization with C 2 puckering deformation. We note that in a real spectrum the low energy window (i.e. below 33000 cm −1 ) will exhibit peaks coming from ESA to singly excited states which remain unresolved with the small active space used in this study. Nevertheless, since these states involve orbitals other than H-1, H, L and L+1, we expect that their signals will blue-shift along the relaxation pathway. We also note that the absolute positions of the fingerprint ESA are likely to change following more accurate computations. Exemplarily, comparison with reference gas-phase calculations indicates that peak 3 appears around 30000 cm −1 . Nevertheless, our findings make peaks 3, 6 and 8 unique for the L a relaxation pathway. In contrast to local excited states, the spectroscopic trace of the CT states is dy- Fig. 4 ESA and SE evolution along the CT relaxation pathway: left) spectral traces at critical points; right) level scheme namic, depending on the energetic fluctuations of the CT energy levels given by the GS dynamics and on the mixing with degenerate transitions. Our calculations indicate that stabilization down to the NUV region alone is not enough to promote a notable fraction of molecules into the CT states upon NUV pumping, as the oscillator strength of CT states remains generally weak, albeit increasing with respect to unstacked conformations. A noteworthy gain is only possible through mixing with the L a band. As a consequence, the CT state would inherit oscillator strength and ESA features of the L a band, whereas the positions of these features would fluctuate with the energy of the CT. In Fig. 4, we focus solely on the spectral characteristics of a pure CT state in order to disentangle in such ideal spectrum its characteristic components that are hard to monitor experimen-tally, due to the weak oscillator strength of pure CT states. In the relaxed GS geometry the lowest CT state (i.e. CT(3'→5')) lies at 6.10 eV with respect to the GS, therefore, pumping occurs with a far-UV (FUV) pulse centered at 49000 cm −1 . Since we are interested in the NUV spectral window, probing is done in the 20000-40000 cm −1 range, where the diagonal GSB of the CT (peak 3) cannot be resolved. An off-diagonal bleach signal at 38900 cm −1 attributed to the L a bleaching (peak 1) is visible. In the FC region the spectral window below the L a bleach is lacking intense ESA, whereas four intense peaks appear at the CT minimum. They originate from strongly red-shifted transitions to mixed doubly excited states which, based on wavefunction analysis, we label CT+L a (3') (peak 4), CT+L a (5') (peak 5), CT+L b (3') (peak 6) and CT+L b (5') (peak 7). Exemplarily, CT+L a (3') denotes a local H→L excitation in the positively charged moiety of the exciplex. The red-shifted broad ESA feature is the fingerprint of the lowest CT states in the aggregate. The strong red-shift indicates that the L a /L b transitions are stabilized along the CT de-excitation pathway in support of the finding that geometrical deformations along these pathways present similarities.

2DNUV spectra for hypothetical deactivation scenarios
At this point we have enough information to generate ideal 2DNUV ultrafast spectra for various scenarios. Thereby, we make use of the experimental observation that in a poly-adenine DNA single strand the decay times (τ 1 = 0.39 ps, τ 2 = 4.3 ps, τ 3 = 182 ps) 70 differ by an order of magnitude. We assume that leaving the FC region, fast and slow decay do not overlap temporally. Thus, in the ps regime only ESA and SE originating from the deactivation channel responsible for the long life-time survive. Fig. 5 shows ultrafast 2D spectra for three hypothetical de-excitation mechanisms. Pumping is done with 733 cm −1 broad pulses (corresponds to a 20 fs Fourier limited pulse) centered at 38000 cm −1 , while 2932 cm −1 broad pulses centered at 34000 cm −1 are used for probing. Plotting is done on a linear scale.

Mechanism A: deactivation exclusively through local excited states.
The first hypothetical mechanism (A, Figure 5 upper panel) that we consider involves the following features: • both exciton states L a (3')+L a (5') and L a (3')-L a (5') are populated upon excitation; • the CT states are not populated in the FC and do not play a role in the de-activation process; • the L a plateau can be reached from each exciton state populated in the FC point through bifurcation of the wavepacket in the femtosecond regime; • due to the plateau on the potential energy surface of L a the wavepacket undergoes spreading in the ps regime ‡ . ‡ The wavepacket spreading along the L a is accounted for by averaging the spectra over 10 snapshots generated by interpolating energies and transition dipole moments between L a (5') 7 and L a (5') 40 , as well as between L a (3') 11 and L a (3') 40 .  According to mechanism A the exciton state L a (3')+L a (5'), which absorbs at 38900 cm −1 dominates the spectrum, the weaker L a (3')-L a (5') state, which absorbs at 38500 cm −1 , introduces asymmetry in the signals. Two ESA bands, one around 33000-35000 cm −1 and another one around 36000-37000 cm −1 are observed. On a femtosecond timescale the SE red-shifts to 31000-32000 cm −1 , the ESA around 33000-35000 cm −1 blue-shifts by 1000 cm −1 and becomes more intense, while the ESA around 36000-37000 cm −1 disappears. On a ps timescale the SE spreads in the Vis part of the spectrum between 26000-32000 cm −1 and consquently its intensity decreases. The ESA around 34000-36000 cm −1 stays compact and retains its intensity. Intuitively, this ESA behavior could be attributed to a compactness of the wavepacket and coherent dynamics in the excited state. This behavior reflects the presense of a flat excited state profile along the L a de-excitation pathway.
3.5.2 Mechanism B: deactivation involving population of CT states through intra-molecular vibrational energy redistribution. The second hypothetical deactivation mechanism (B, Figure 5 middle panel) involves the following features: • both exciton states L a (3')+L a (5') and L a (3')-L a (5') are populated upon excitation; • the CT states are not populated in the FC; • the L a plateau can be reached from each exciton state populated in the FC point through bifurcation of the wavepacket in the femtosecond regime; • when reaching the L a plateau of each base the wavepacket bifurcates again; most of the L a population decays rapidly (i.e. in the few-ps regime) to the GS, while a non negligible fraction is transfered to a CT channel through intra-molecular vibrational energy redistribution where it gets trapped in the corresponding minimum; • as a result of the observations of Conti et al. CT(3'→5') is accessible only from the L a (5') plateau, while CT(5'→3') is accessible only from the L a (3') plateau.
The femtosecond photophysics of the first and second mechanisms are identical. With the majority of L a population decaying back to the GS in the ps regime the overall spectral intensity decreases. The decay of the L a population is coded in the partial recovery of the GSB. The population transfer to the CT shows in the broad ESA band between 31000 cm −1 and 37000 cm −1 which emerges along the L a (3')-L a (5') and L a (3')+L a (5') traces. Computations by Conti et al. suggest that mechanisms A and B are correlated, whereby slowing down the relaxation along the flat L a plateau (mechanism A) facilitates the back-and-forward population transfer between the L a and CT channels. Provided this mechanism is confirmed experimentally, it is intriguing whether it is a coherent process showing through oscillations between the L a and CT traces.
1-18 | 13 3.5.3 Mechanism C: deactivation involving population of CT states in the FC region. The third hypothetical deactivation mechanism (C, Figure 5 lower panel) involves the following features: • in 50% of the snapshots the CT(3'→5') is stabilized dynamically below the L a region to ∼38000 cm −1 and can be directly populated; § note that the percentage and the energetic postion of the CT(3'→5') state are arbitrary despite some theoretical support 64 ; • the CT(3'→5') state mixes with the L a band in the FC region, borrows oscillator strength and adopts its ESA features ¶ ; • local transitions (L a , D * 1 (3'/5'), D * 2 (3'/5') and D * 3 (5')) are not affected by the deformations (shortening of the inter-base distance) which stabilize the CT state; • the L a plateau can be reached from each exciton state populated in the FC point through bifurcation of the wavepacket in the femtosecond regime; • the L a population decays rapidly (i.e. in the few-ps regime) to the GS, the CT population is trapped in the local minimum on the excited state.
Although this mechanism is not supported by QM/MM computations 31,32,58 , other groups have considered the possibility to populate low lying CT states directly in the FC region 8,64 . In this regard we would like to draw attention to a recent paper by Rohlfing and co-workers which argues that aqueous solvation may lower CT states by more than 1 eV based on many-body Green's function computations 64 . We believe that only a detailed analysis of the conformational dynamics considering thermal fluctuations instead of analyzing a single representative relaxed geometry will give a definite answer to this pending question. At the FC point the CT trace appears in the long-wavelength region of the NUV absorption band. Due to mixing with the L a band L a features appear as weak ESA along the CT trace at ∼ 35500 cm −1 and ∼ 38000 cm −1 , blue-shifted by 1000 cm −1 with respect to the intense ESA along the L a (3')+L a (5') trace. In the femtosecond timescale the L a channel relaxes to its plateau, while the CT channel relaxes to its local minimum Min CT(3 →5 ) and the CT fingerprints are clearly recognized. In the ps regime only the CT fingerprints survive. The fs population recovery is again coded in the GSB of both the L a and CT bands. Note, that neither GSB recovers completely as there is still population in the excited state. Comparison between the second and thrid mechanisms demonstrates the sensitivity of 2D spectroscopy in resolving population transfer (B). Furthermore, 1D pump-probe spectroscopy, in which the 2D spectra colapse along ω 1 , will hardly discriminate the ps signatures of these two mechanisms. § This is accomplished by red-shifting the energy of the CT(3'→5') state, as well as the energies of the mixed doubly excited states CT+L a (3'), CT+L a (5'), CT+L b (3') and CT+L b (5') (giving rise to peaks 4-7 in Fig. 4) by 8000 cm −1 . ¶ This is accomplished by assigning finite transition dipole moments from the CT(3'→5') state to the local doubly excited states accessible from the L a states (D * 1 (3'/5'), D * 2 (3'/5') and D * 3 (5') in Figure  3).
In this contribution we provided a framework for simulating time-resolved 2D electronic spectroscopy experiments relying on accurate ab initio methods. The SOS//QM/MM protocol requires transition energies of singly and doubly excited states together with the transition dipole moments among them. While dynamics simulations are necessary to cover all feasible de-excitation pathways and to collect enough information for simulating non-homogeneous broadenings, already static computations (i.e. geometry optimizations, minimum energy path calculations) may provide information on the temporal evolution of individual deactivation channels and may be used to generate ideal 2D spectra in support of hypothetical scenarios, derived from the mechanistic picture. To assess the validity of these ideal spectra a thorough modeling on top of mixed quantum-classical dynamics simulations is required. The SOS//QM/MM protocol is a hybrid approach and as such it can profit from the progress in three different fields of research: QM/MM dynamics, excited state computations and non-linear spectroscopy. Potential improvements range from employing polarizable force fields 71 or orbital-free embedding 72 to treat the environment in the QM/MM dynamics via using more flexible wavefunction-based techniques for computing excited states (e.g. generalized active space SCF 73 , the density matrix renormalization group approach 74 , the n-electron valence state perturbation theory 75 ) to using real pulse shapes and accounting for coherence evolution when computing the nonlinear response of the system 28 . The adenine-adenine stacked dimer embedded in a explicitly solvated doublestranded DNA model was chosen to demonstrate the application of the protocol in practice. Characteristic points along the first excited state potential energy surface were used to generate time-dependent traces of the local L a and CT excimer channels. Subsequently, the 2DNUV spectra were simulated for three deexcitation mechanisms: a deactivation exclusively along the L a surface (mechanism A), a deactivation involving population of CT states through intra-molecular vibrational energy redistribution (mechanism B) and a deactivation associated with populating both the L a and the CT states in the FC region. A more detailed study will follow using more characteristic points along the potential energy surface, applying more accurate computational methods (e.g. increasing the active space), considering the role of the L b band, including coherence evolution and parametrized inhomogeneous broadenings, as well as applying a kinetic model to simulate population transfer along and between deactivation channels. Fig. S1. List of the valence π-orbitals of adenine together with different single electron transitions which determine the low energy part of the absorption spectrum. Nearly degenerate transitions that are likely to mix are combined in groups denoted by a color code.

Reference calculation for adenine monomer in gas-phase
A reference calculation of adenine monomer in the gas-phase at a state-of-the-art SA-CASSCF/SS-CASPT2 level was performed to assess the validity of the approximations used in the main article. The protocol for the reference calculations has been first reported in ref. 1. The gas-phase geometry was optimized at single state CASSCF level with a full π-valence active space of twelve electrons in ten orbitals (i.e. CAS (12,10)) under CS-symmetry with the generally contracted ANO-L basis set 2 , whereas the following contraction scheme is adopted: C,O/[4s3p2d] and H/[2s1p]. It is well know that excited states described with the general ANO basis set are contaminated by spurious Rydberg character that can affect excited state energies and properties like transition dipole moments 3,4 . To account for this the basis sets was augmented with a set of 8s, 8p and 8d uncontracted Rydberg-type basis functions positioned at the center of nuclear charge of the lowest two cationic states 5,6 in all excited state calculations. We point out that we intentionally used a set of uncontracted Rydberg-type basis functions instead of the usually utilized 1s, 1p, 1d contracted Rydberg-type set 7 in order to properly describe Rydberg states up to 10.50 eV. The orbital exponents for the uncontracted Rydberg basis functions were obtained from ref. 8. The uneven description of electronic correlation for valence and Rydberg states at CASSCF level is the cause for valence-Rydberg mixing 7 which affects electronic properties of valence states. Well described Rydberg orbitals do not interact with the valence orbitals 5,9 , therefore, the following scheme was applied to resolve valence-Rydberg orbital mixing: An ten electrons / nine orbitals active space (i.e. CAS (10,9)) comprising all valence orbitals besides the completely bonding orbital was constructed. The active space was augmented with a Rydberg orbital of A''-symmetry (the Rydberg orbitals of A' symmetry are not considered further), thus forming CAS(10,10). Next a state-averaged calculation was performed where it was paid attention that at least one of the optimized states had Rydberg character. Subsequently, the active space orbitals were localized using the Pipek procedure 10 in order to prevent spurious valence-Rydberg orbital mixing. Finally, the perfectly formed Rydberg orbital was deleted from the molecular orbital list. This procedure was performed 24 times until all Rydberg orbitals were removed and, consequently, no Rydberg states were found among the first 25 excited states. Once the valence orbitals were clean from Rydberg contaminations the lowest occupied orbital was added to the active space, thus, forming the full π-valence active space of twelve π-electrons in ten orbitals (i.e. CAS (12,10)) and the first 25 states with A'-symmetry were computed. To test the stability of the calculation the active space was systematically increased by four, eight and twelve additional extravalence secondary orbitals. The idea is that if the contribution of extravalence orbitals to the overall wavefunctions of the described states is negligible the CASPT2 excitation energies should not differ if the extravalence orbitals are treated variationally or perturbatively. If this is not true, than the orbitals must be included in the AS. A convergence was reached with twelve secondy orbtials, resulting in the active space RAS(0,0|12,10|2,12), where each of the three restricted active subspaces are characterized by the maximum number of simultaneously excited electrons and the number of orbitals. To reduce the computational effort only configurations with up to two excited electrons in these virtual orbitals were added to the list of configuration state functions. The contributions of the remaining configurations were treated perturbationally with the multiconfigurational counterpart of the Møller-Plesset method denoted as CASPT2/RASPT2 11 . All virtual orbitals were correlated in the perturbation procedure except for the 24 removed Rydberg orbitals with A''-symmetry. An IPEA shift of 0.0 a.u. 12 and an imaginary shift of 0.2 a.u. 13 were used. Transition dipole moments were calculated at SA-CASSCF/RASSCF level. The RICD approximation was used to speed up the calculation of two-electron integrals 14 .   20A' (D + )