Solvent effect on the energetics of proton coupled electron transfer in guanine-cytosine pair in chloroform by mixed explicit and implicit solvation models

Accepted Manuscript Solvent effect on the energetics of Proton Coupled Electron Transfer in GuanineCytosine pair in chloroform by mixed explicit and implicit solvation models Lara Martinez-Fernandez, Giacomo Prampolini, Javier Cerezo, Yanli Liu, Fabrizio Santoro, Roberto Improta PII: DOI: Reference:

S0301-0104(18)30454-3 https://doi.org/10.1016/j.chemphys.2018.07.012 CHEMPH 10076

To appear in:

Chemical Physics

Received Date: Revised Date: Accepted Date:

2 May 2018 26 June 2018 17 July 2018

Please cite this article as: L. Martinez-Fernandez, G. Prampolini, J. Cerezo, Y. Liu, F. Santoro, R. Improta, Solvent effect on the energetics of Proton Coupled Electron Transfer in Guanine-Cytosine pair in chloroform by mixed explicit and implicit solvation models, Chemical Physics (2018), doi: https://doi.org/10.1016/j.chemphys. 2018.07.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Solvent effect on the energetics of Proton Coupled Electron Transfer in Guanine-Cytosine pair in chloroform by mixed explicit and implicit solvation models. Lara Martinez-Fernandez,a Giacomo Prampolini,b Javier Cerezo ,c Yanli Liu,d Fabrizio Santoro,b and Roberto Improta*ae a) LIDYL, CEA, CNRS, Université Paris-Saclay, France b) CNR−Consiglio Nazionale delle Ricerche, Istituto di Chimica dei Composti Organo Metallici (ICCOM-CNR), SS di Pisa, Area della Ricerca, via G. Moruzzi 1, I-56124 Pisa, Italy c) Departamento de Química Física, Universidad de Murcia, 30100 Murcia, Spain d) School of Physics and Optoelectronics Engineering, Ludong University, 264025 Yantai, People's Republic of China e) Consiglio Nazionale delle Ricerche, Istituto di Biostrutture e Bioimmagini, 80136 Naples, Italy * Corresponding author E-mail: [email protected] Abstract The Watson-Crick hydrogen bonded pair formed by deoxyGuanosine and deoxyCytidine (GC) in chloroform has been analysed by classical Molecular Dynamics simulations, which shows the existence of several fairly stable solute/solvent hydrogen bonds. Time-Dependent Density Functional Theory (TD-DFT) calculations with M052X functional, including solvation effect by a mixed continuum/discrete model, show that solvent stabilizes G C Charge Transfer excited states (CTGC). The Proton Transfer (PT) processes that can occur in CTGC have been mapped by TD-DFT combining State Specific and Linear Response implementations of the Polarizable Continuum Model. For the first time, the effect of explicit solute/solvent interactions on the PT is considered by studying models containing up to three CHCl3 molecules. Our study shows that PT from the N1 atom of G to the N3 of C is exoergonic also in solution but, at variance with what observed in gas-phase, a stable minimum is predicted for CTGC state in chloroform. 1. Introduction Domcke and collaborators in a series of seminal papers, have shown that Proton Coupled Electron Transfer (PCET) can provide a very effective excited state deactivation channel for hydrogen bonded (H-bonded) pairs formed, in the gas phase, by DNA bases or even by some simple base models.[1-6] In the guaninecytosine (GC) pair, connected by Watson-Crick (WC) Hydrogen Bonds (HB), the lowest energy bright excited state can give access to an inter-molecular charge transfer state (CTGC).[3,4] CTGC leads to a barrierless proton transfer (PT) reaction, which yields a biradical species (G-H) (C+H) (Scheme 1). Such a process is more efficient in the WC form than in other alternative H-bonded arrangements, since the relative stability of CTGC is larger for WC pairs.[3] By invoking this PCET process, which in the following we shall label as “Domcke mechanism”, Domcke and Sobolewski provided a convincing explanation of the much shorter excited state lifetime of WC GC pair, witnessed by its very broad absorption spectrum.[7] A similar

1

mechanism has been shown to be operative in aminopurine hydrogen bonded dimers,[5,6] and it could be possible in the dimers formed by adenine (A) thymine (T)[8] and by anionic[9,10] oxoguanine (with A and C).[11] Later computational studies basically confirmed the results of Domcke and Sobolewski, showing that in the gas phase PCET processes can be operative between H-bonded nucleobases dimers, giving account of the experimental ultrashort excited state lifetime.[12-15] The picture is instead more controversial for what concerns the results obtained in the condensed phase and, even more, within DNA duplexes. In a moderately polar solvent as chloroform, by using bulky substituents on the sugar, it is possible to induce the formation of WC pairs containing GC and AT nucleosides, giving the possibility to verify if PCET processes occur also in solution.[16-19] For what concerns the chloroform solution containing G and C, Schwalb and Temps found that the fluorescence (FU) lifetime (<0.3 ps) in GC WC H-bonded pairs (and in GG pairs) was much shorter than in the monomer, for excitation wavelength comprised between 296 (close to the band origin) and 262 nm.[19] They ascribed this result to the involvement of a PCET process as that depicted in Scheme 1. A similar conclusion has been recently reached by analysing Transient Absorption (TA) and Transient Infra-Red (TIR) spectra of a GC WC pair[16] thanks to the individuation of the spectral signature of (G-H) radicals. Based on such identification, it has been proposed that the (G-H)(C+H) diradical species is formed in ≤ 40 fs while a lifetime of 2.9 ps was assigned to the decay of (G-H)(C+H) to the ground electronic state. Bleach recovery signal was observed in 7.2 ps, due to the necessary full vibrational cooling. According to Röttger at al.[16] the sub-ps decay found in FU experiments is thus related to the formation of the GC CT state and not to ground state recovery. Furthermore, the TIR spectra suggested also the formation (ca 10% of the photoexcited population) of a stable photoproduct, produced by double PT (PT2GC Scheme 1). When exciting at 266 nm, the PCET route was estimated to have a quantum yield ≥ 0.6, which is, however, significantly smaller at 290 nm. This difference could explain why, following excitation at 283 nm, the broad TA spectra of GC pairs dimers (in 270-700 nm range) is similar to that of the separated monomers, suggesting a minor involvement of PCET in the ultrafast dynamics.[20] No spectral signature of PCET was instead found for AT WC pairs in chloroform.[17]

2

GC

PT2GC

PTGC Ÿ

CTGC δ+

Ÿ

δ-

PTGCNH

2

Ÿ

Ÿ

Scheme 1: General scheme of the possible charge transfer and proton transfer mechanisms triggered by UV absorption in GC pair. Please note that PTGC and PTGCNH2 could in principle be formed by H atom transfer from the other bright excited states.  For what concerns the duplex, many excited state decay routes are possible, making much more difficult unravelling their different spectral signatures. In particular, intra-strand CT excited states are also possible and, according to the available computational studies, they are more stable than the inter-strand CT states, i.e. those occurring between WC pairs, which the present study focuses on.[21,22] Actually, there are now many converging indications, both experimental and computational, that a significant part of the excited state population can decay to intra-strand CT states[23-34] that can successively trigger inter-strand PT processes.[23,24] On the other hand, although computational studies suggest that Domcke mechanism is possible also within a duplex,[35] no unambiguous experimental proof exists that it is operative in DNA.[36] However, recent TR studies on DNA miniduplexes show that, when stacking is prevented, an additional decay route, characterized by ground state recovery in 4.6 ps,[37] is present and it has been tentatively ascribed to inter-strand PCET. 3

It is clear, then, that, whereas the involvement of Domcke process in the gas phase can be considered well assessed, the question is still open in solution or in DNA. In any case, experiments clearly show that even a moderately polar solvent as chloroform dramatically affects the PCET reaction, increasing the excited state lifetime and decreasing the quantum yield. The PT reaction, leading to a neutral diradical species, quenches the strong dipole moment of the G+C- ion pair, making a significant dependence on solvent likely, both from the static and the dynamic points of view. The significant changes of the electron density associated to the PCET reaction make, indeed, an accurate treatment of dynamical solvation effects fundamental. Full equilibration of solvent degrees of freedom to the changes in the electron density caused by the electronic excitation requires, in fact, a finite time[38-42] (on the ps time scale), which depends on the solvent. Actually, several computational studies, exploiting the State-Specific (SS) Polarizable Continuum Model (PCM) in combination with TD-DFT (Time-Dependent Density Functional Theory),[43,44] indicate that inclusion of the solvent significantly changes the energetic of the PCET reaction with respect to the gas phase, even in a moderately polar solvent as chloroform.[20,45] When solvent is considered at the equilibrium with the excited state density (equilibrium calculations) (i) the driving force towards the PCET reaction decreases, (ii) an energy barrier appears on the associated path and (iii) the energy difference between the di-radical minimum and the Conical intersection (CoIn) with the ground state increases.[21] On the one hand, comparison of non-equilibrium and equilibrium SS-PCM/TD-DFT predictions usually provides a reliable description of dynamical solvent effects on transitions with a significant amount of CT character,[43,44] and the indications of the calculations are consistent with an increase of the time associated to ground state recovery via the PCET route. On the other hand, it is clear that strong and specific solute-solvent interactions make desirable the adoption of an explicit solvation model, taking into account the molecular nature of the solvent. As a first step towards a fully dynamical study of PCET in the condensed phase, we here report a detailed characterization of the energetics of the PT processes occurring in the lowest-energy inter-strand CT state in a GC WC pair in chloroform solution. We combine Molecular Dynamics (MD) classical simulations and PCM/(TD)-DFT calculations, using M052X as reference density functional. After having characterized the solvation shell of GC in chloroform by MD simulations, we shall analyse how the solvent affects the stability of the different excited states in the Franck Condon (FC) region, focussing on the role played by explicit solute/solvent interactions. By combining SS-PCM and LRPCM-TD-DFT calculations we shall map the different PT processes that can occur in the CTGC, showing that a stable minimum is found for this state when the computational model includes explicitly three chloroform molecules (in the following indicated as CHCl3) or dynamical solvation effects are included at the SSPCM/TD-DFT level.

4

2.Computational details 2.1 MD simulations All simulations were performed with the Gromacs 5.1.1 package,[46,47] in the NPT ensemble. Unless otherwise stated, temperature (T, 298 K) and pressure (P, 1 atm) were kept constant through the v-rescale and Parrinello-Rahaman schemes,[48,49] employing τT = 0.1 ps and τP = 5 ps, as coupling constants. For short-range inter-molecular interactions, a cut-off radius of 14 Å was employed, while for long-range electrostatic interactions, the Particle Mesh Ewald (PME) scheme was adopted. In all runs a 0.1 fs time step was used, since no constraint was applied on the fast-vibrating stretching coordinates.

C6 C5 N1

C4 N8

N3

H8

O11

C2 H1 N1

O7 H10

C6 C5 C8 N10 C4 N3 N9 C2

S0 minimum

PTGC-min

N7

CTGC-min

PTGCNH2-min

Figure 1: Schematic drawing and atom labelling of the most important structures discussed in this study. See text and Scheme 1 for more details on the different species. The starting geometry of the GC pair, shown in Figure 1, was obtained through QM calculations, carried out at Density Functional Theory (DFT) level, using the M052X functional and the 6-31+G(d,p) basis set. In all MD simulations, to describe both GC intra- and inter-molecular interactions, the GAFF force field along with AM1-BCC charges [50] was adopted, and the topology prepared through the acpype parser of the antechamber tool.[51] Chloroform was described through the GAFF force-field. All parameters and the

5

starting structure, containing 1000 CHCl3 molecules, were downloaded from the virtualchem.org website (http://virtualchemistry.org).[52] A preliminary equilibration of the bulk liquid phase was carried out for a total of 10 ns, at 298 K and 1 atm. The GC pair DFT optimized geometry was solvated into the chloroform liquid through the genbox tool of the GROMACS package.[46,47] The resulting configuration, composed by one GC pair and 989 CHCl3 molecules, was employed as starting arrangement for a series of equilibration runs, carried out as outlined in the following: 

A short equilibration run (200 ps) to observe the first impact of the solvent and check for any artifact due to the initial solvation procedure.



A long equilibration run (5 ns) to let the solvent equilibrate correctly around the solute. To allow the density to reach equilibrium more quickly, the Berendsen scheme[53] was used for both thermostat and barostat.



A long production run (8 ns) to retrieve the snapshots within the correct ensemble (NPT) distribution. To this end, the v-rescale and Parrinello-Rahaman thermostat/barostat combination[48,49] was preferred. During the production run, the MD trajectory of the whole system was saved every 10 ps.

The solvation structure and the specific solvent-solute interactions were investigated over the stored NPT trajectories by computing the Spatial Distribution Functions (SDF) and the atom-atom pair-correlation functions, g(r), defined as:

where  is a target site of the GC pair and  an atom of the kth CHCl3 molecule, while the mean value stands for a double average over all frames and solvent molecules. 2.2. QM calculations DFT has been used to optimize the ground electronic state (S0) of the GC pair, and TD-DFT to characterize the excited states. M052X[54,55] has been our reference functional, because it provides good performance when treating dispersions interactions (which are important to describe the clusters containing CHCl3 molecules) and a fairly accurate description of CT transitions,[54,55]. Due to these properties, M052X has been already successfully applied to the study of DNA photophysics and photochemistry.[21,56-62] The picture of the PCET process provided by M052X has been double-checked by using two additional 6

functionals: the long-range corrected CAM-B3LYP[63] and M062X[64]. Our previous studies [??] shows that in the gas phase our TD-DFT calculations provide a picture of the PCET reaction very similar to that obtained by WF-based methods, supporting their reliability for describing this process. 6-31G+(d,p) basis set has been used in geometry optimizations and in excited calculations, complemented by some test calculations by using the 6-311+G(2d,2p) basis set. The excited state calculations on the 100 structures issuing from MD simulations have been performed on the G and C nucleosides by using the smaller 6-31G(d) basis set. Bulk solvent effects have been included by PCM,[39,65] performing geometry optimizations by its linear response (LR) version, for which analytical gradients are available.[66] The energetics of the PCET process have been refined by resorting to SS implementation (SS-PCM) TD-DFT calculations,[43,44] which provide a much more reliable description of dynamical solvation effects, especially when dealing with CT transitions. In our implementation, starting from LR TD-DFT calculation the state specific reaction field is computed using the electron density of the state of interest. In the following step, a TD DFT calculation is performed in the presence of this first set of polarization charges, providing an updated excited state density and, consequently, a new set of polarization charges. This iterative procedure is continued until convergence on the reaction field is achieved. Both non-equilibrium (neq) and equilibrium (eq) time-regimes have been considered. In the former, slow solvation degrees of freedom remains in equilibrium with the ground electronic density, whereas only the electronic polarization (fast solvation degrees of freedom) is equilibrium with the excited-state electron density. In the latter time-regime, solvent is fully equilibrated with the excited state electronic density. The PCET process has been mapped also by some test calculations by using the Corrected Linear response method.[67] All the calculations have been performed by the Gaussian09 Program.[68] 3. Results 3.1. MD simulations The starting structure of the CG pair, shown in Figure 1, indicates an almost coplanar disposition of the aromatic backbones of the two moieties, which are held together by a network of three HB, established between the O11-H8, H1-N3 and H10-O7 pairs. This disposition is maintained upon solvation: from a first visual inspection of the MD production run at 298 K and 1 atm, it clearly appears that the solvent embeds the 7

whole pair, which conserves its planarity while the HB network is only seldom broken and always quickly reestablished (see Supporting Information for a more detailed analysis). MD simulations also provide insights into the microscopic structure of the first solvation shell, useful to assess which regions of the GC pair are more prone to interact with the solvent. In Figure 2, the SDF computed over the MD production runs are reported in three different views and shows that the highest solvent density is in the vicinity of the HB network that holds together the GC pair.

Figure 2: SDFs for the GC pair in CHCl3 solvent. Explicit ball&stick representation was adopted for all atoms of the solute, whereas the density of H and Cl atoms of the solvent is displayed in white and green, respectively. a) side view; b) top view; c) end view.

Figure 3: Top: g(r) computed between the O11 (left), O7 (center) and N3 (right) GC pair atoms and the H (blue) or Cl (green) of any CHCl3 molecule. The end of the first neighbor peak is evidenced with a vertical red dotted line. Bottom: number of H or Cl atoms of the solvent within the first neighbor peak of atom  (O11, O7 or N3), obtained by proper integration of the pair correlation function. GC atom labels refer to Figure 1. 8

This preliminary analysis calls for further investigation on the solvation structure, focusing in particular on the region with a higher solvent density, which also coincides with the area more involved in the PCET. To this end, the atomic pair correlation functions, g(r), have been computed between the most relevant GC sites and the solvent atoms. Inspection of Figure 3, where all the resulting curves are displayed, shows two general trends. First, for all considered interactions sites, it is clear from the position of the first peak that the solvent molecules approach the solute by pointing their H atoms toward the bases, rather than Cl ones. The most pronounced first neighbor peaks involve the external oxygen atoms, indicating that C-O7 and GO11 are by far the GC sites most exposed to the solvent. Hydrogen atoms of the neighboring CHCl3 molecules can be found within a 2 to 2.5 Å radius from each Oxygen, suggesting stable solute-solvent HB. Another HB interaction, though weaker than those with the Oxygen atom, involves the Lone Pair (LP) of the C-N3 atom. The g(r) pair correlation functions can be further exploited to extract information on the average number of solvent molecules, which can be found at short interaction distances from the G-O11, C-O7 and C-N3 sites. Indeed, the first peak maximum is located in both cases around 2.3 Å, a distance compatible with an additional HB-like interaction, established between the GC pair and one or more neighboring solvent molecules. Such HB could in principle strongly affect the PCET process, and will therefore be taken explicitly into account in QM calculations. An estimate of the number of solvent molecules strongly interacting with the GC pair can be extracted by integrating the g(r) functions, as shown in the bottom panels of Figure 3. It appears that one to two molecules are able to come in close contact with the G-O11 atom, whereas one CHCl3 molecule in average can be found within a radius of 3 Å from C-O7. In summary, therefore, two or three CHCl3 molecules should be explicitly included in the QM calculations, to take into account more accurately the possible perturbations induced by the solvent on the PCET process. Finally, despite the inner GC HB seems to be more screened from specific interactions with the surrounding solvent, a chloroform molecule, H-bonded to C-N3, an additional weak HB might be established between a chloroform molecule and C-N3. 3.2. The excited electronic states in the FC region Based on the outcome of our MD simulations, we have then analysed the excited state behaviour of GC pairs in the FC region, by using two complementary approaches. In the first approach, we have studied 9methylguanine-1methylcytosine Watson-Crick dimer (in the following denoted simply as GC), including a different number of explicit CHCl3 molecules (from 0 to 3), whose position has been chosen on the ground of our MD analysis. In the largest model two CHCl3 have been placed so to form a HB with the two Carbonyl Groups of C (Chl ) and G (Chl ), whereas the third can interact with the Nitrogen LP of C-N3 (Chl , see Figure 4 and S6). This latter molecule is not considered in the model including only two CHCl3 molecules. Finally, we have investigated the model where only Chl  is explicitly considered. For all the models two 9

slightly different local minima (I and II) have been investigated (see also SI), and re-optimized in chloroform at the PCM/M052X/6-31+G(d,p) level of theory. Table 1 reports the Vertical Absorption Energy obtained at the PCM/TD-M052X/6-31+G(d,p)//PCM/M052X/6-31+G(d,p) level. For isolated C in chloroform only one absorption band is recorded in experiment below 250 nm, and its maximum falls at 280 nm.[18] According to our calculations this peak is due to a HOMO-LUMO * excitation, similar to that already described in the gas phase and in other solvents. Two features are found in the experimental absorption spectrum of G falling at 275 and 260 nm and both more intense than the peak of C. For the isolated monomers in chloroform a pure continuum model of the solvent can capture the main experimental features. Indeed, our calculations predict only one absorption peak for C, less intense and slightly red-shifted with respect to the two transitions of G, which corresponds to the La and Lb * transitions.[69] According to our calculations, formation of the GC dimer leads to significant changes of the monomer electronic transitions. For instance, some mixing between the two lowest energy electronic transitions of G and C is observed and the lowest energy excited state S1 is mainly localized on G, whereas S2 is mainly on C, reversing the relative energy of the two transitions in the monomers (Table 1). This result can be explained by the sensitivity of the C-1* to strong HB, which is well known from experimental and computational studies in H-bonded solvents.[58,70] The effect of WC pairing on the relative energy of G-La and C-1* transition is confirmed by CC2 calculations[71] and CASPT2 calculations[12] in the gas phase. A third intense transition can be associated to the Lb transition of G (hereafter G-Lb) whereas, the lowest energy G C CT state (hereafter CTGC), that can be described as an excitation from the G HOMO to the C LUMO, corresponds to S4 and is ca 0.15 blue-shifted with respect to S3.

Figure 4: Structures and selected solute-solvent distance in the S0 minimum (a) and in the excited state GCCT-min (b), according to LR-PCM/TD-M052X/6-31+G(d,p) geometry optimizations.

10

From the quantitative point of view the computed band are uniformly (0.7 eV) blue-shift with respect to the experimental one. Besides to the inaccuracy of M052X, this discrepancy is due in part (roughly 0.1~0.2 eV, based on our previous experience [72]) to the absence of vibronic and thermal effects in our calculations and part to the size of the basis set and to the absence of explicit solvent molecules (see below), accounting roughly for an additional 0.1~0.2 eV shift. The PCM only picture is basically confirmed when CHCl3 molecules are explicitly included in the calculation. The most significant difference with respect the ‘PCM only’ results concerns a slight stabilization of the CTGC excited state, which can become more stable than G-Lb. Actually, there is one model where these two states are extremely mixed and share their oscillator strength (column 5 in Table 1). The stabilization of CTGC is also related to a small mixing with close-lying solute-solvent CT state involving G and one CHCl3 molecule, which is obviously absent in the ‘PCM only’ calculations. Table 1: Vertical Absorption Energies (in eV) and Oscillator Strength (in parentheses) of the Four (Five) Lowest Excited States Computed at the LR-PCM/TD-M052X/6-31+G(d,p) level for GC in Chloroform solution. Two different clusters (I and II) are shown for 3 CHCl3 (see SI). PCM only

PCM+

PCM+

PCM+

PCM+

1 CHCl3 I

2 CHCl3 I

3 CHCl3 I

3 CHCl3 II

5.06(0.12)

5.08(0.09)

5.02(0.13)

4.96(0.13)

4.95(0.12)

5.22(0.28)

5.23(0.27)

5.20(0.25)

5.15(0.22)

5.14(0.21)

5.51(0.58)

5.50(0.53)

5.50(0.51)

5.55 (0.31)

5.66(0.01)

5.59(0.03)

5.65(0.02)

5.70(0.01)

5.73(0.00)

Description

G- La

C- * 1

a

5.50(0.53)

5.46 (0.23)

a

5.44(0.02)

5.73 (0.00)

5.71(0.00)

G- Lb

CTGC G CHCl3 CT

1

Notes: Isolated G (PCM only) G-La 5.21 (0.21), G-Lb 5.58 (0.49); Isolated C (PCM only) C- * 5.14 (0.19); a) strongly mixed excited states.

In the second approach, we have computed at the full QM level the six lowest excited states on 100 snapshots issuing from MD simulations, explicitly including all the CHCl3 molecules (between 13 and 19) within 3 Å with respect to any of the GC atoms and taking into account bulk solvent effect at the PCM level (see Computational details). The simulated absorption spectrum (obtained broadening each transition with a Gaussian with half width half maximum= 0.05eV) is shown in Figure 5. The most intense band falls at 5.15 11

eV/240 nm and is preceded by a small shoulder in the 250-260 nm region. S1, S2, and S3 in a lesser extent, are responsible for this absorption band. A second less intense peak is present below 230 nm, which shows contributions from all the remaining excited states (S3, S4, S5 and S6). A more detailed analysis shows that S1 has mainly G-La character (in agreement with the results obtained on QM optimized geometries shown in Table 1) and falls within the 4.90-5.10 eV range for 60% of the snapshots. However, for some snapshots this state is S2 being less stable than the one with C-1* character. The assignment of the other bright state, G-Lb, is more complex since for most of the geometries it is strongly coupled with GCHCl3 CT excited states. Indeed, its oscillator strength varies from 0.1 to 0.5 and is not always the third excited state (S3) but, depending on the degree of mixing, it can fall higher in energy (up to S6). The same conclusions hold for the CTGC state, i.e. we found different state ordering and energies whereas its oscillator strength is, in general, rather low (<0.1). Furthermore, as described above for the GCHCl3, the CTGC appears also often coupled with the bright states within an energy difference lower than 0.4 eV. 1

Norm. Abs.

0.8 0.6 0.4 0.2 0 200

210

220

230

240

250

260

270

wavelength / nm Figure 5: Simulated absorption spectra at the PCM/TD-M052X/6-31G(d) level of theory using the 100 snapshots extracted from the MD (‘MD/QM’) in solid black line, broadening each transition with a Gaussian function with hwhm=0.05 eV. In red (pure QM) the simulated spectra for the QM clusters PCM (solid) and PCM+3 CHCl3 I (dashed), with the same level of theory and hwhm=0.30 eV. The comparison of the ‘pure QM’ and ‘ MD/QM’ (i.e. arising from the snapshot calculations) spectra reveals two main differences (Figure 5). First, there is a significant energy shift of the absorption maxima, redshifted by ~ 0.6 eV in the MD/QM spectrum. Second, in the latter spectrum the relative intensity of the band falling at 225nm, which should be related to the G-Lb, strongly decreases with respect to the ‘pure QM’ spectrum. These results can be in principle due to a combination of different factors: (i) the use of the 12

GAFF FF, instead of the DFT functionals, to describe the geometry of the solute (ii) the explicit inclusion of many solvent molecules in the MD/QM calculations, (iii) the inclusion in the latter of thermal effects. In order to get some additional insights on this issue, we have computed the spectrum of 20 snapshots issuing from MD simulations at the ‘PCM only’ level, i.e. without including any explicit solvent molecule. The comparison made in Figure S4 shows that the low energy bands are not significantly affected by inclusion of CHCl3 molecules, but for a very weak red-shift of the absorption maximum, which is also slightly broader. This result suggests that the use of MM geometries could lead to an artificial red-shift of the absorption maximum,

confirming

the

results

obtained

when

computing

the

absorption

spectrum

of

5methylCytosine.[70] The underestimation of the G-Lb intensity could also be due to some inaccuracies of the MM geometries. Actually, the GC pair geometries issuing from a simple MM geometry optimization performed with the GAFF parameters and those averaged over the MD simulations are rather different from those optimized at the PCM/M052X level (see Figure S5), both concerning the monomer internal structures and the HB network interconnecting the two bases. As far as the internal geometries are concerned, it might be worth noticing that negligible differences arise between MM optimized and MD averages structures, but significant distortions are found with respect to the DFT optimized arrangement. In particular, some C-C double bonds are significantly longer in the MM and MD structures, leading to a smaller energy gap between occupied and virtual orbitals, and, therefore, to red-shifted spectra. The larger differences stand however in the inter-base HB, which are found significantly larger in the average geometry arising from simulation than in the two optimized structures (DFT and MM). In this case, along with the effect of the solvent, thermal effects could also play a role. To further investigate this issue, we performed a supplementary MD run, on the isolated GC pair in the NVT ensemble at 298 K. The comparison between the inter-pair HB distances (see Table S2), as obtained according to DFT and MM optimization or MD simulations, performed in vacuo or in chloroform solvent is indeed instructive, revealing that thermal effect rather than the presence of the solvent are the main responsible for the elongation of the HB with respect to the DFT optimized structure. On the other hand, the blue-wing of the spectrum is significantly affected by inclusion of CHCl3 molecules; being the ‘PCM only’ spectra broader and more intense. This comparison confirms the large coupling between * and CT states involving solvent molecules for the majority of snapshots, which could be mirrored in this lower intensity. Finally, several trials increasing the number of excited states until 10 showed that no additional state with significant oscillator strength affects the investigated spectral region.

3.3. Mapping the PCET process in continuum solvation model 13

Excited state geometry optimizations of the CTGC for GC, after a crossing to the close-lying bright excited states, leads directly, without any energy barrier, to the PT reaction between G-N1 and C-N3, converging the di-radical minimum shown in Figure 1 (hereafter PTGC-min). Constrained geometry optimizations, preventing the transfer of G-H1 to C-N3 and of the G-NH10 amino proton to C-O7, lead to the CTGC pseudominimum (CTGC-min) where G adopts the geometry typical of a G cation and C that of C anion, which is characterized, inter-alia, by strong pyramidalization of C its amino moiety. This feature probably explains why in CTGC-min the two bases undergo a noticeable rotation around the inter-molecular axis, reaching a twisted arrangement, which allows preserving the HB. The HB lengths G-H1/C-N3 and G-NH10/C-O7 are very short, confirming the tendency of G+ to undergo PT towards C-. Actually, if the C-NH10 distance is not blocked, geometry optimization transfers this proton to the C-C2-O7 carbonyl group, finding another diradical minimum (hereafter PTGCNH2-min), which is a stationary point of the S1 Potential Energy Surface (PES). As CTGC-min, PTGC-min and PTGCNH2-min exhibit a twisted geometry. CTGC-min is less stable than PTGCmin and PTGCNH2 by 0.66 eV and by 0.17 eV, respectively. These energy differences are not affected by an increase of the basis set, being 0.66 eV and 0.17 eV also at the PCM-TD-MO52X/6-311+G(2d,2p) level of theory. Interestingly, in the gas phase, the relative stability of CTGC-min decreases by 0.2 eV so that it is less stable than PTGC-min and PTGCNH2 by 0.9 eV and 0.39 eV, respectively. As discussed more in detail below, CTGC-min is more polar than PTGC-min and PTGCNH2, explaining why is relatively stabilized also by a moderately polar solvent as chloroform. The most significant difference between the results obtained in the gas phase and in chloroform concerns the energy gap from S0. In chloroform, at the LR-PCM/TD-M052X level, PTGC-min is 2.4 eV less stable than S0, whereas in the gas phase the energy difference is smaller than 0.4 eV, suggesting the presence of an easily accessible CoIn with the ground state, fully in line with the conclusions of previous computational studies in the gas phase.[4] Including solvent effect on the PCET reaction, thus, decreases the accessibility of the decay channel to S0 from PTGC-min. A relaxed scan in PCM (LR calculations at the neq level) performed by increasing the G-N1-H1 and G-N10H10 distances (see Figure 6) provides indications coherent with what got from full geometry optimizations: starting from CTGC-min, both PT reactions are barrierless and exo-ergonic. Detailed inspection of Figure 6 shows that a small energy barrier (≤ 2 kcal/mol) separates PTGC-min from PTGCNH2-min and that there is a wide region of the PES, energetically accessible, where both G N-H distances are elongated, suggesting the possibility of a double PT reaction. Actually, the structures corresponding to a Double Proton Transfer (hereafter GC-DPT, upper right corner in Figure 6; notice that more precisely the transfer involves a proton and a hydrogen atom) are very close in energy to PTGCNH2-min, though no clear minimum is present for this structure in the excited states. This picture is fully consistent with the thorough analysis performed by Sauri 14

et al. at the CASSCF/CASPT2 level in the gas phase,[12] considering that in the gas phase PTGC-min and PTGCNH2-min minima are much closer to the CoIn with S0 than in chloroform solution. As a consequence, the transfer of a hydrogen atom from PTGC-min or PTGCNH2-min leads to the crossing to S0 and, actually, GC-DPT corresponds to a stable ground state minimum.

Figure 6: Contour plot (energy in kcal/mol) obtained from a relaxed scan along the G-N1-H1 and G-N10-H10 distances. PCM/TD-M052X/6-31+G(d,p) calculations on GC dimer in chloroform. Please note that PTGC-min would be further stabilized by an increase of the G-N1-H1 distance. As shown also in previous papers,[43,44,73] LR-PCM/TD-DFT is not suitable to describe electronic transitions involving a large variation of the electron density and underestimates solvent effect on strongly ‘dipolar’ excited states as GCCT. As a consequence, we have recomputed the energy of the different points of the scan at the State Specific PCM/TD-M052X level, both in the non-equilibrium and in the equilibrium (see Computational details) time-regimes (see Figure 7), focussing on the transfer of the N1-H proton from G to C, which is the most favoured according to our LR-PCM/TD-M052X analysis. As shown in Figure 7, already at the SS-PCM/TD-M052X neq level, where only the fast solvation degrees of freedom are in equilibrium with the excited state density, CTGC-min (involving a transfer of 0.75 a.u. electron charge from G to C) is stabilized with respect to S0 by 0.7 eV; it also become closer in energy to PTGC-min, whose dipole moment is much smaller (7.0 Debye versus 13.6 Debye). The energy difference between CTGC-min and PTGC-min indeed decreases down to 0.14 eV (neq level) and 0.15 eV (eq level).

15

Furthermore a small energy barrier (0.15 eV) is associated to the PCET process. At the eq level the energy barrier becomes twice as large (0.30 eV) and both CTGC-min and PTGC-min are significantly closer to S0. It should be noted that, when the solvent is in equilibrium with S1, S0 is destabilized, since only fast solvent degrees of freedom are equilibrated with its own electron density. PTGC-min remains more stable than CTGCmin, though by only 0.15 eV.

Figure 7: Evolution of the CTGC excited state energy as a function of the G(N1-H1) bond distance, computed in chloroform by LR-PCM/TD-M052X/6-31+G(d,p) calculations and SS-PCM/TD-M052X/6-31+G(d,p). In the latter case the curves obtained within the non-equilibrium (neq) and equilibrium (eq) time-regime are shown. As a next step of our analysis we have verified how the energetics of PCET is affected by (i) the choice of the density functional and (ii) the radii used to build the cavity in PCM calculations. As shown in the Supporting Information, CAM-B3LYP and M062X provide a picture very close to that of M052X. For what concerns the PCM parameters, we performed some test calculations by changing the atomic radii from UFF to UAKS, which lead to a smaller solute cavity, its volume being 30% smaller than when using UFF radii. For UAKS calculations, CTGC-min and PTGC-min are almost isoenergetic (energy difference < 0.01 eV), in line with a larger stabilization of the more dipolar form in the presence of a smaller cavity, i.e. when the solvent is ‘closer’ to the solute. The picture predicted by SS-PCM/TD-M052X calculations is qualitatively confirmed by test calculations adopting the corrected linear response model (see SI), though no energy barrier is predicted for the PCET at this level of theory. 16

3.4 Describing PCET by a Mixed continuum-discrete solvation model In this section we analyse how the inclusion of explicit chloroform molecules affects the PCET process. To this aim we have optimized CTGC in the presence of 3 CHCl3 molecules (see Figure 4). At difference with what happens in ‘PCM only’ calculations, CTGC-min is predicted to be a stationary point of the PES. The two molecular planes exhibit a significant ‘bending’, allowing the maximization of the interaction with the CHCl3 molecules. These latter are remarkably displaced with respect to the coordination geometry adopted in the S0 minimum, allowing the maximization of the HB strength between CHCl3 (acting as HB donor) and the C base, which is negatively charged in CTGC-min. For example, the CHCl3 molecule denoted as  in Figure 4, which is H-bonded with the G-C6-O11 group in the S0 minimum, shifts towards C-N3 at the CTGC-min. On the same time, CHCl3  approaches the Lone Pair of the C-NH2 amino group. Interestingly, CHCl3 , which is already H-bonded to the C-C2-O7 carbonyl group, tilts towards the G+ moiety in order to allow a stronger interaction with the chlorine atoms. Geometry optimization of PTGCmin in the presence of 3 CHCl3 molecules indicates that the solvation shell is more similar to that found in the S0 minimum (see SI) with respect to that optimized for CTGC-min. PTGC-min is confirmed to be more stable than CTGC-min, but the energy gap, computed at the LR-PCM/TD-M052X level, decreases by ca 0.1 eV up to 0.55 eV. 4. Concluding Remarks In this paper we have analysed how chloroform affects the energetic of the PCET process involving the inter-base CT state in a GC WC dimer. First of all, we have performed an extensive analysis of the dynamical behaviour of the GC dimer in chloroform solution, by exploiting classical MD. Our simulations give interesting insights on the solvation shell of GC and provide a more solid basis for the subsequent full QM analysis. At least one CHCl3 molecule is predicted to form a stable HB with G-O10 atom of the G carbonyl group, and structures with two H-bonded molecules to this atom are also common. Analogously, one stable solute/solvent HB involves C-O7, while also C-N3 frequently acts as HB acceptor with a solvent molecule. These results provide very useful guidelines for building smaller computational models to be studied by full QM geometry optimizations both in the ground and in the CT excited state. As detailed below, we have indeed examined different computational models, including up to 3 CHCl3 molecules of the first solvation shell. For what concerns the analysis of the FC region, our calculations highlights that, independently of the nature of the embedding medium, the relative energy of the lowest energy transitions of G and C is reversed upon WC pairing. While for the isolated molecules the lowest energy absorption maximum of C is red-shifted with respect to that of G, in GC the lowest energy excited state is mainly localized on G. This result should be considered when analysing the dependence of the excited state dynamics on the 17

excitation wavelength.[16] It indicates in fact that it cannot be taken for granted that decreasing the excitation wavelength from 260 nm to 290 nm, i.e. on the red-wing of the absorption band, C is preferentially photoexcited, as it would happen for the isolated bases. On the other hand, our calculations confirm that at 260 nm, i.e. on the blue-wing of the absorption band, G (due to the Lb contribution) absorbs much more than C. Actually, for what concerns the dependence of the excited state dynamics on the excitation wavelength, a critical issue is the relative position of the GCCT state with respect to the bright states. Our calculations indicate that GCCT falls in the proximity of G-Lb and that the inclusion of solvent effect can lead to significant mixing of these two states. The largest PCET yield observed at 260 nm could thus be related to an increase of the population of the GCCT exited state. Both our ‘pure QM’ and ‘MD/QM’ approaches indicate that the explicit inclusion of CHCl3 molecules in the calculations does not have a significant effect on the red-wing of the computed absorption spectrum but can affect its blue-wing. In particular, the GCCT state is predicted to be stabilized and more mixed with the G-Lb bright state, making more likely a population transfer. Calculations including explicitly CHCl3 solvent molecules also suggest that the possible role of solute/solvent CT states should be investigated more in detail, since GCHCl3 CT states falls close in energy to the other excited states in the FC region. Experiments do not show any significant involvement of solvent molecules in the excited state dynamics, but, on the other hand, GCHCl3 CT could lead to the formation of very reactive chlorine radicals, which could interfere in the excited state reactivity. Once the system decays to GCCT our PCM/TD-M052X calculations predict that G-N1-HC-N3 PT from the G+ to the C- base is exoergonic and that also in solution this PT is favoured over G-N10-H10C-O7 PT reaction. This picture does not depend on the functional and is analogous to what obtained in the gas phase using different computational methods. On the other hand, our analysis confirms that solvent strongly affects the PCET reaction, especially when solvation models, as State Specific PCM, are used. First of all, in solution the minima resulting from PT transfer, both PTGC-min and PTGCNH2-min, are much more distant in energy from the CoIn with S0 than in the gas phase. Furthermore, the relative stability of CTGC-min increases, since, as already anticipated in the introduction, an increase of the polarity of the embedding medium stabilizes the ionic G+ and C- species. This effect is operative already at the LR-PCM/TD-M052X level, although this method underestimates the stability of CT minima. When adopting SS-PCM/TD-DFT calculations, indeed, the relative stability of CTGC-min and PTGC-min becomes much more similar and a small energy barrier is associated to the PT process. This prediction is confirmed also by the calculations studying the effect of explicit CHCl3 molecules on the PCET energetics in GC. In the GC+3 CHCl3 model, CTGC-min is predicted to be a stable minimum of the PES already at the LR-PCM/TD-M052X level. The solvation shell is sensitive to the presence of the G+ and C- ionic species and rearranges according to the ‘new’ electronic distribution of the CT state. In particular, HB between the solvent molecules and the C- moiety, where CHCl3 mainly acts as HB 18

donor to the Lone Pair of the N and O atoms of the base, are strengthened. From a different perspective, HB to the N3 atom (partially negatively charged in the CT minimum) decreases its basicity, i.e. the tendency of acquiring an extra proton from the G+ moiety. This result provides a caveat for the description of this PT process by using a continuum solvation model and motivates future works explicitly including a larger number of solvent molecules in the description of PCET PES. . Additionally, it cannot be taken for granted that using the same solvent cavity for the different excited states does not affect the description of this excited state process. This conclusion is confirmed by our test calculations using a smaller PCM cavity, which strongly stabilizes CTGC-min. PTGC-min is, however, confirmed to be more stable than CTGC-min also when three CHCl3 molecules are included in the calculation. This result suggests that the occurrence of the PCET process in solution is likely, especially for excitations in the proximity of G-Lb state, when the possibility of populating the CT state is higher, and this prediction is in line with the experimental indications. On the other hand, our calculations suggest that the characteristic time associated to the PCET process (< 40 fs) on the ground of the experimental results is too short. It is likely that, at least for certain solvent configurations, part of the photoexcited population remains trapped for sometime in the CT minimum. Actually, as also acknowledged by Rottger at al., the absorption spectrum of the G+ cation and of the G-H radical (i.e. the product of the PT reaction) are very similar.[74,75] it is thus possible that both CT and PT states contribute to the spectral features appearing on the 50 fs time-scale. The presence of an energy barrier on the PT reaction path, together with the energy gap between PTGC-min and S0 can give instead account of the ‘longer’ excited state lifetime (2.9 ps) assigned to the CT/PT state. In this respect it is noteworthy that the energy gap between PTGC-min and S0 strongly decreases for eq calculations; the 2.9 ps time-constant is thus consistent with our suggestion that a full equilibration of solvation degrees of freedom to the excited state density, which occurs on a few ps time scale, is crucial to reach the CoIn with S0. Furthermore, a small energy barrier can contribute to the larger PT yield observed for shorter excitation wavelengths, i.e. when more energy is deposited on the system. Though providing new insights on solvent effect on the Domcke mechanism, our approach overlooks the first step of this reaction, i.e. the population transfer from the bright states towards CGCT. In other words, we have implicitly assumed a stepwise process, where the systems first decay to the CT states and, then, PT is triggered. Actually, it would be also possible than the diradical species PTGC-min is directly reached by the bright excited states via the transfer of a hydrogen atom.[12] Taking this first step of the reaction into the proper account would also make more difficult to describe dynamical solvation effects. In the simplest picture, soon after photoabsorption the slow solvation degrees of freedom are still in equilibrium with S0, whereas the fast ones with the populated bright excited state (we can label generically as SB). Population transfer to the CT, occurring while the slow solvent degrees of freedom are equilibrating with SB, would 19

make the situation more complex, even before the PT reaction leads to a dramatic change of the electron density. It is thus therefore clear that only a real a dynamical treatment, suitably including the coupling between intra-molecular and solvation degrees of freedom, can provide a solid interpretation of the experimental results. This paper constitutes the first, preliminary, step in this direction, which will be the goal of our future studies in the field. Similar considerations apply also to the study of PCET within DNA. It is clear that environmental effects would play a critical role in the photoactivated reactions, considering as ’environment’ not only the solvent but the entire phospho-deoxyribose backbone, whose conformational fluctuations is expected to strongly modulate the associated PES. On the other hand, the strong anisotropic nature of the duplexes, together with the necessity of suitably describing structural rearrangements occurring on multiple time-scales, probably make unavoidable resorting to purely atomistic computational approaches. Acknowledgments RI thanks Program D’Alembert of the University Paris-Saclay for financial support. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

B. Marchetti, T.N.V. Karsili, M.N.R. Ashfold, W. Domcke, Phys. Chem. Chem. Phys. 18 (2016) 20007. A.L. Sobolewski, W. Domcke, J. Phys. Chem. A 111 (2007) 11725. A.L. Sobolewski, W. Domcke, C. Hättig, Proc. Natl. Acad. Sci. 102 (2005) 17903. A.L. Sobolewski, W. Domcke, Phys. Chem. Chem. Phys. 6 (2004) 2763. T. Schultz, E. Samoylova, W. Radloff, I.V. Hertel, A.L. Sobolewski, W. Domcke, Science 306 (2004) 1765. A.L. Sobolewski, W. Domcke, Chem. Phys. 294 (2003) 73. A. Abo-Riziq, L. Grace, E. Nir, M. Kabelac, P. Hobza, M.S. de Vries, Proc. Natl. Acad. Sci. 102 (2005) 20. S. Perun, A.L. Sobolewski, W. Domcke, J. Phys. Chem. A 110 (2006) 9031. A.-O. Colson, B. Besler, D.M. Close, M.D. Sevilla, J. Phys. Chem. 96 (1992) 661. A. Kumar, M.D. Sevilla, Chem. Rev. 110 (2010) 7002. X. Wu, N.T. Karsili, W. Domcke, Molecules 22 (2017) 135. V. Sauri, J.P. Gobbo, J.J. Serrano-Pérez, M. Lundberg, P.B. Coto, L. Serrano-Andrés, A.C. Borin, R. Lindh, M. Merchán, D. Roca-Sanjuán, J. Chem. Theor. Comput. 9 (2013) 481. J.P. Gobbo, V. Saurí, D. Roca-Sanjuán, L. Serrano-Andrés, M. Merchán, A.C. Borin, J. Phys. Chem. B 116 (2012) 4089. P.R.L. Markwick, N.L. Doltsinis, J. Schlitter, J. Chem. Phys. 126 (2007) 045104. M.K. Shukla, J. Leszczynski, J. Phys. Chem. A 106 (2002) 4709. K. Röttger, H.J.B. Marroux, M.P. Grubb, P.M. Coulter, H. Böhnke, A.S. Henderson, M.C. Galan, F. Temps, A.J. Orr-Ewing, G.M. Roberts, Ang. Chem. Int. Ed. 54 (2015) 14719. K. Röttger, H.J.B. Marroux, A.F.M. Chemin, E. Elsdon, T.A.A. Oliver, S.T.G. Street, A.S. Henderson, M.C. Galan, A.J. Orr-Ewing, G.M. Roberts, J. Phys. Chem. B 121 (2017) 4448. N.K. Schwalb, F. Temps, J. Am. Chem. Soc. 129 (2007) 9272. N.K. Schwalb, T. Michalak, F. Temps, J. Phys. Chem. B 113 (2009) 16365. L. Biemann, S.A. Kovalenko, K. Kleinermanns, R. Mahrwald, M. Markert, R. Improta, J. Am. Chem. Soc. 133 (2011) 19664. L. Martinez-Fernandez, R. Improta, Faraday Discuss. 207 (2018) 199. 20

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

C. Ko, S. Hammes-Schiffer, J. Phys. Chem. Lett. 4 (2013) 2540. Y. Zhang, K. de La Harpe, A.A. Beckstead, L. Martínez-Fernández, R. Improta, B. Kohler, J. Phys. Chem. Lett. 7 (2016) 950. Y. Zhang, K. de La Harpe, A.A. Beckstead, R. Improta, B. Kohler, J. Am. Chem. Soc. 137 (2015) 7059. P.M. Keane, M. Wojdyla, G.W. Doorley, J.M. Kelly, A.W. Parker, I.P. Clark, G.M. Greetham, M. Towrie, L.M. Magno, S.J. Quinn, Chem. Commun. 50 (2014) 2990. G.W. Doorley, M. Wojdyla, G.W. Watson, M. Towrie, A.W. Parker, J.M. Kelly, S.J. Quinn, J. Phys. Chem. Lett. 4 (2013) 2739. M.C. Stuhldreier, F. Temps, Faraday Discuss. 163 (2013) 173. D.B. Bucher, B.M. Pilles, T. Carell, W. Zinth, Proc. Natl. Acad. Sci. 111 (2014) 4369. B.M. Pilles, D.B. Bucher, L.Z. Liu, P. Gilch, W. Zinth, W.J. Schreier, Chem. Commun. 50 (2014) 15623. Y. Zhang, J. Dood, A.A. Beckstead, X.B. Li, K.V. Nguyen, C.J. Burrows, R. Improta, B. Kohler, Proc. Natl. Acad. Sci. 111 (2014) 11612. Y. Zhang, J. Dood, A.A. Beckstead, X.B. Li, K.V. Nguyen, C.J. Burrows, R. Improta, B. Kohler, J. Phys. Chem. B 119 (2015) 7491. G.W. Doorley, D.A. McGovern, M.W. George, M. Towrie, A.W. Parker, J.M. Kelly, S.J. Quinn, Ang. Chem. Int. Ed. 48 (2009) 123. D.B. Bucher, A. Schlueter, T. Carell, W. Zinth, Ang. Chem. Int. Ed. 53 (2014) 11366. L. Martinez-Fernandez, Y. Zhang, K. de La Harpe, A.A. Beckstead, B. Kohler, R. Improta, Phys. Chem. Chem. Phys 18 (2016) 21241. G. Groenhof, L.V. Schäfer, M. Boggio-Pasqua, M. Goette, H. Grubmüller, M.A. Robb, J. Am. Chem. Soc. 129 (2007) 6812. F.-A. Miannay, Á. Bányász, T. Gustavsson, D. Markovitsi, J. Am. Chem. Soc. 129 (2007) 14574. Y. Zhang, X.-B. Li, A.M. Fleming, J. Dood, A.A. Beckstead, A.M. Orendt, C.J. Burrows, B. Kohler, J. Am. Chem. Soc. 138 (2016) 7395. R. Improta, V. Barone, F. Santoro, J. Phys. Chem. B 111 (2007) 14080. J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev. 105 (2005) 2999. M. Maroncelli, G.R. Fleming, J. Chem. Phys. 86 (1987) 6221. M.L. Horng, J.A. Gardecki, A. Papazyan, M. Maroncelli, The Journal of Physical Chemistry 99 (1995) 17311. E. Pluhařová, P. Slavíček, P. Jungwirth, Acc. Chem. Res. 48 (2015) 1209. R. Improta, V. Barone, G. Scalmani, M.J. Frisch, J. Chem. Phys. 125 (2006) 054103. R. Improta, G. Scalmani, M.J. Frisch, V. Barone, J. Chem. Phys. 127 (2007) 074504. M. Dargiewicz, M. Biczysko, R. Improta, V. Barone, Phys. Chem. Chem. Phys. 14 (2012) 8981. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, E. Mark Alan, J.C. Berendsen Herman, J. Comput. Chem. 26 (2005) 1701. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M.R. Shirts, J.C. Smith, P.M. Kasson, D. van der Spoel, B. Hess, E. Lindahl, Bioinformatics 29 (2013) 845. M. Parrinello, A. Rahman, Journal of Applied Physics 52 (1981) 7182. G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 126 (2007) 014101. J. Wang, M. Wolf Romain, W. Caldwell James, A. Kollman Peter, A. Case David, J. Comput. Chem. 25 (2004) 1157. A.W. Sousa da Silva, W.F. Vranken, BMC Research Notes 5 (2012) 367. C. Caleman, P.J. van Maaren, M. Hong, J.S. Hub, L.T. Costa, D. van der Spoel, J. Chem. Theor. Comput. 8 (2012) 61. H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. Y. Zhao, N.E. Schultz, D.G. Truhlar, J. Chem. Theor. Comput. 2 (2006) 364. Y. Zhao, D.G. Truhlar, Acc. Chem. Res. 41 (2008) 157. I. Conti, L. Martinez Fernandez, L. Esposito, S. Hofinger, A. Nenov, m. Garavelli, R. Improta, Chem. Eur. J. (2017) n/a. L. Martinez-Fernandez, R. Improta, Photochem. Photobiol. Sci. 16 (2017) 1277. 21

[58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75]

L. Martinez-Fernandez, A.J. Pepino, J. Segarra-Martí, J. Jovaisaite, I. Vayá, A. Nenov, D. Markovitsi, T. Gustavsson, A. Banyasz, M. Garavelli, R. Improta, J. Am. Chem. Soc. 139 (2017) 7780. R. Improta, J. Phys. Chem. B 116 (2012) 14261. A. Banyasz, L. Esposito, T. Douki, M. Perron, R. Improta, D. Markovitsi, J. Phys. Chem. B 120 (2016) 4232. A. Banyasz, L. Martinez-Fernandez, T. Ketola, A. Muñoz-Losa, L. Esposito, D. Markovitsi, R. Improta, J. Phys. Chem. Lett. 7 (2016) 2020. L. Esposito, A. Banyasz, T. Douki, M. Perron, D. Markovitsi, R. Improta, J. Am. Chem. Soc. 136 (2014) 10838. T. Yanai, D.P. Tew, N.C. Handy, Chemical Physics Letters 393 (2004) 51. Y. Zhao, D.G. Truhlar, Theor. Chem. Acc. 120 (2008) 215. S. Miertuš, E. Scrocco, J. Tomasi, Chem. Phys. 55 (1981) 117. G. Scalmani, M.J. Frisch, B. Mennucci, J. Tomasi, R. Cammi, V. Barone, Journal of Chemical Physics 124 (2006) 094107. M. Caricato, B. Mennucci, J. Tomasi, J. Chem. Phys. 124 (2006) 124520. M.J. Frisch et al., ( Gaussian 09, revision A.1; Gaussian, Inc.,Wallingford, CT, 2009). R. Improta, F. Santoro, L. Blancafort, Chem. Rev. 116 (2016) 3540. L. Martínez-Fernández, A.J. Pepino, J. Segarra-Martí, A. Banyasz, M. Garavelli, R. Improta, J. Chem. Theor. Comput. 12 (2016) 4430. S. Yamazaki, T. Taketsugu, Phys. Chem. Chem. Phys. 14 (2012) 8866. F.J.A. Ferrer, J. Cerezo, E. Stendardo, R. Improta, F. Santoro, J. Chem. Theor. Comput. 9 (2013) 2072. R. Improta, Phys. Chem. Chem. Phys. 10 (2008) 2656. A. Banyasz, L. Martínez-Fernández, C. Balty, M. Perron, T. Douki, R. Improta, D. Markovitsi, J. Am. Chem. Soc. 139 (2017) 10561. L.P. Candeias, S. Steenken, J. Am. Chem. Soc. 111 (1989) 1094.

22