Photochemistry of Thymine in Protic Polar Nanomeric Droplets Using Electrostatic Embeding TD-DFT/MM

Thymine photochemistry is important for understanding DNA photodamage. In the gas phase, thymine undergoes a fast non-radiative decay from S 2 to S 1 . In the S 1 state, it gets trapped for several picoseconds until returning to the ground-state S 0 . Here, we explore the electrostatic effects of nanomeric droplets of methanol and water on the excited states of thymine. For this purpose, we develop and implement an electrostatic embedding TD-DFT/MM method based on a QM/MM coupling defined through electrostatic potential fitting charges. We show that both in methanol and water, the mechanism is similar to the gas phase. The solvent molecules participate in defining the branching plane of S 0 /S 1 intersection and have a negligible effect on the S 1 /S 2 intersection. Despite the wrong topology of the ground/excited state intersections, electrostatic embedding TD-DFT/MM allows for a fast exploration of the potential energy surfaces and a qualitative picture of the photophysics of thymine in solvent droplets.


Introduction
Photodamage of deoxyribonucleic acid (DNA) originating from ultraviolet (UV) sunlight radiation is one of the main sources of skin cancer [1]. UV light can induce irreparable damage to the genetic code, despite the many natural protection and correction mechanisms that exist [2]. The most basic protection against UV radiation is found at the nucleobase level. Indeed, it has been shown that isolated guanine, adenine, thymine and cytosine undergo ultrafast non-radiative decays back to the ground state [3][4][5][6][7]. Photoexcitation with UV light brings the nucleobase in a singlet excited state of ππ * character, which relaxes to lower-energy dark states of nπ * character (n is a lone-pair electron of a heteroatom in the system) and finally back to the ground state.
The nucleobases decay mechanisms depend strongly on the embedding environment [8][9][10][11][12]. Indeed, when nucleobases are stacked in a DNA polymer, internal non-radiative decay competes with several charge-and proton transfers which ultimately lead to cyclobutane dimers or 6-4 photoadducts [13][14][15][16][17][18], at the origin of DNA mutations. Electronic structure combined with quantum dynamics simulations are especially suited to determine the mechanisms of deactivation in the excited state. However, there is a lack of methodologies that can treat nucleobases in their native environment at a reasonable computational cost. Especial care has to be taken at the level of theory taken for the electronic structure since the energetic order and the critical points of the photoreaction highly depend on the amount of electron correlation included in the simulations [19]. For example, three different photochemical mechanisms have been inferred from theoretical simulations of thymine: S 1 trapping (fast S 2 and slow S 1 decays) [20], S 2 trapping (slow S 2 and fast S 1 decays) [21], and S 1 /S 2 trapping (slow S 1 and S 2 decays) [22].
Recently, the effect of the electronic structure on the photochemistry of thymine in the gas phase has been studied, showing that thymine undergoes a fast non-radiative decay from S 2 to S 1 state thanks to a barrierless access to a conical intersection followed by an S 1 -trapping, which appears because of a barrier to reach the conical intersection back to the ground state [19]. Here, we extend the study of Ref. [19] to determine the effect of protic solvent nanomeric droplets on the photochemistry of thymine. For this purpose, we develop and implement a new quantum mechanics/molecular mechanics (QM/MM) hybrid electrostatic embedding method combining time-dependent density-functional theory (TD-DFT) [23], and classical force fields via an electrostatic potential fitting (ESPF) interaction Hamiltonian [24]. The ESPF method has been recently implemented for the analytic energy, gradient and Hessian for DFT including charge conservation and grid derivatives [25][26][27]. This method has been shown to successfully represented the electrostatic interaction between complex electrostatic environments such as proteins or solvents [28,29]. The new implementation extends these developments to the excited state, formulating the energy and analytic gradient of TD-DFT/MM required for performing both static (critical point search) and dynamic (semi-classical quantum dynamics) photochemical studies.
We apply the TD-DFT/MM method implemented here to the photochemistry of thymine in water and methanol nanomeric droplets. Despite being only qualitatively correct, TD-DFT has been shown to successfully describe the potential energy surfaces of thymine in the aqueous solution described via continuum models [9]. TD-DFT allows for a fast exploration of potential energy surface (PES) for which excited state mechanisms can be easily extracted. Indeed, we can search for minimum energy structure and approximate conical intersections via the branching update method [30]. Conical intersections are the critical points at the origin of most ultrafast non-radiative decay mechanisms, alike to the transition states in the ground state reactivity. Optimization of conical intersections is thus fundamental for determining the static picture photochemical mechanisms.
This study is organized as follows: in Section 2 the DFT/MM and TD-DFT/MM methods are described, in Section 3 the implementation and computational details are detailed, and in Section 4 the results on the photochemistry and conical intersection topology of thymine in the gas phase and in methanol and water droplets are shown. The conclusions are detailed in Section 5

Methodology
The basics of the ESPF QM/MM electrostatic embedding coupling with ground-state self-consistent field single-determinant methods have been published and discussed in previous literature [24,26,27]. Here we followed the notation and method described in Ref. [27], in which a major reformulation of ESPF was established. Here, a formulation is presented based on ESPF that extends the QM/MM coupling to excited states. Without loss of generality, we kept the derivations using the Tamm-Dancoff approximation of the linear response TD-DFT equations (i.e., decoupled excitation and de-excitation parts of the linear response function) [31].

Electrostatic Embedding Coupling in ESPF DFT/MM
In the electrostatic potential fitting (ESPF) method, the QM and MM subsystems are put in contact through a purely electrostatic interaction, that can be written as (1) In this equation, φ A is the MM electrostatic potential on the QM center A, where q FF i is the MM partial charge defined from the force field, and R A and R i are the Cartesian coordinates of QM and MM atoms respectively. Hereafter, we will consider only classical force fields, that is, q FF i has a fixed value without polarization. The q A in Equation (1) is the partial charge on the QM atom A, which is defined as where Z A is the atomic charge of atom A, N AO is the number of atomic orbitals, P µν is the one-electron atomic-orbital density matrix, and Q A,µν is the electron charge population operator on QM center A. This operator is obtained by a fitting procedure, in which the electrostatic integrals on a grid are fitted into operators on the QM centers. This accounts to solving the equation Using the short-hand notation T k,A = |r k − R A | −1 and V k,µν = µ| 1 |r−r k | |ν = dr |r−r k | , we can formally write the electron charge operators as Deriving Equation (1) with respect to the density matrix, we arrive at the ESPF Hamiltonian, which is given simply by Due to the finite size of the grid defined in Equation (4), the charge conservation condition, is not exactly satisfied. In Ref. [27], I showed that this can be imposed at the level of operators by deriving the charge conservation condition with respect to the density matrix, in which S µν = µ|ν is the atomic-orbital overlap matrix. This charge conservation in operator form can be imposed in the Hamiltonian directly, leading to where we introduced the average external potential, The Hamiltonian defined in Equation (9) can thus be added to the gas phase QM Fock operator (F 0 ) to define an electrostatic embedding formulation. The total energy of the system then becomes in which E MM is the MM energy contribution and Z A is the atomic number of atom A. The gradient of Equation (11) has two blocks with a different expression, depending on the type of atom (QM or MM) with respect to the energy is derived, in which we symbolized the derivatives E x = ∂E/∂x, in which x andx represent a Cartesian coordinate of a QM or MM atom respectively. In the QM block, we substituted the derivative of the atomic-orbital density matrix with the trace of the derivative of the overlap and the energy-weighted matrix W = PFP † . Since the ESPF charge operator only depends on the coordinates of the grid points and the QM atoms, the derivative of the ESPF operator with respect to MM atom perturbations is simply given by in which we defined Q A = Tr AO [PQ A ]. Thus, the gradient becomes

Electrostatic Embedding Coupling in ESPF TDA-TDDFT/MM
In TDA-TDDFT/MM, the excited state energy of state I is obtained by adding the excitation energy ω I to Equation (11). The excitation energy is obtained by solving the eigenvalue problem, , in which f Hxc = δv H /δρ + δv xc /δρ is the Hartree-exchange-correlation kernel obtained by deriving the corresponding potentials with respect to the electronic density [32]. The energy gradient is efficiently obtained by the Lagrangian approach by Furche and Ahlrichs [33]. In this approach, TDA-TDDFT/MM excitation energy gradient for state I is written as In this equation, P I and Γ I are the one-and two-electron relaxed transition densities to state I, W I = P I FP † I and q A,I = Z A − Tr[P I Q A ] (for the definition of relaxed transition densities, see Ref. [33]). The relaxed densities are obtained after solving a set of coupled-perturbed equations (the Z-vector equations), which only contain two-electron terms and thus not affected by the ESPF one-electron interaction Hamiltonian [33,34].

Branching Update Method
In the branching update method [30], the optimization is performed by minimizing the average energy E = E I + E J /2 + E J − E I 2 using the gradient in which the average gradient is given by ∇E av = ∇E I + ∇E J /2 and the gradient difference is given by ∇E di f f = (∇E I − ∇E 2 )/2. The projector T projects out the branching plane vectors from the average gradient, and is defined as, in which ∇E orth is an vector orthogonal to ∇E di f f . The tilde indicates that these vectors are normalized. The ESPF QM/MM, difference gradient vector between two excited states is defined as in which we defined ∆Q A,I J = Q A,J − Q A,I . The particular case of the excited/ground state gradient, the expression for the gradient difference vector is simply the excitation energy of excited state I As for the orthogonal gradient, it is estimated by a simple update procedure, in which k indicates the optimization step.

Computational Details
The methods described here were implemented in a development version of GAMESS-US (14 FEB 2018 R1) interfaced with a modified version of Tinker 6.3.3 [35,36]. For the computation of ESPF atomic charges, an atom-centred Lebedev grid of 50 angular points and nine radial increments on each atom. Points inside the Van der Waals radii were excluded from the grid. The one-electron ESPF Hamiltonian and the derivative were added accordingly to the one-electron Hamiltonian and derivative parts of the one-electron Hamiltonian of the Fock operator. The gradients of the ground and excited states on MM atoms were computed using the corresponding ESPF atomic charges computed respectively for the density matrix and the relaxed density matrix. The environment was either optimized using the quadratic convergence optimization method of GAMESS-US or using microiterations. All calculations were performed at the TD-BH&HLYP/6-31G* using the Tamm-Dancoff approximation for the thymine atoms [37][38][39], and Amber94 for the water or methanol atoms [40] .

Results and Discussion
In this section, the photochemistry of thymine was studied in gas phase and in nanometic droplets of water and methanol. The purpose was to illustrate the potential of the methodology described in Section 2 rather than performing a full mechanistic study, that would require performing quantum dynamics. The nanomeric droplets were created using the PACKMOL package [41], and subsequently optimized to the minimum energy structure of the internal QM/MM energy. For methanol and water droplets, 100 and 300 molecules respectively were included in the simulations.

Photochemistry of Thymine in Gas Phase
The potential energy surface representing the key points of the photochemistry from the lowest singlet energy states of thymine in the gas phase is summarized in Figure 1. The first singlet (S 1 ) state was described mainly as a single excitation from a lone-pair orbital of oxygen (labeled n, with contributions from oxygen 4a mainly and oxygen 2a to a lesser extent) to a π * orbital delocalized mainly over carbons 4, 5 and 6. The S 2 state was described essentially as a π to π * transition. The photochemical pathways of deactivation of thymine have been discussed in depth in the literature [5]. Here, we focused on the topology of the potential energy surfaces (PES) and the conical intersection seams found at the TD-BH&HLYP/6-31G* level. The PES in the figure was obtained by a geodesic interpolation between the optimized structures of three conical intersections and two state minima [42]. The geodesic path represents a minimum energy path connecting the different critical points of the photoreaction, for which we can extract an idea of the photochemical mechanism. From the S 0 state, thymine was excited to the bright S 2 state of ππ * nature. There was no stable minimum for the S 2 state predicted at the TD-DFT level. Rather, the S 2 minimum was found at the seam of conical intersection with S 1 , of nπ * nature. From the S 1 /S 2 intersection, a branching occurs. On the left side of the graph, the path towards the intersection between the ground state and the ππ * state is shown. This intersection could only be reached by surmounting a barrier of ca. 0.13 eV, and indeed it has been shown as a minor path of deactivation of thymine in the gas phase [5]. Most of the wavepacket would rather follow a minimum energy path towards the nπ * state minimum. This minimum was long-lived, since the corresponding conical intersection was 0.8 eV higher in energy. The intersection geometries, vectors and PES around the intersection point for intersections S 0 /nπ * , S 0 /ππ * and nπ * /ππ * are represented in Figure 2. As seen from the PES along with the non-adiabatic coupling, both intersections involving the S 0 state were sloped conical intersections, while the nπ * /ππ * was peaked. In the S 0 /nπ * intersection, the C 4 -O 4a bond got perpendicular to the thymine plane, while in the S 0 /nπ * intersection the C 5 -C 5a bond was perpendicular to the thymine plane. The nπ * /ππ * was a boat-like structure, in which the bonds N 3 -H and C 6 -H were parallel but displaced with respect to the plane of thymine. The wrong intersection topology for the S 0 /nπ * state was observed in the cut along the gradient difference vector, as noted previously for TDA-TD-DFT in the gas phase, in which a negative excitation energy was leading to the "true" ground state [43]. Indeed, TDA-TD-DFT led to diabatic surfaces rather than adiabatic states along the gradient difference vector.

Photochemistry of Thymine in Methanol Droplet
The photochemistry from the lowest energy singlet states of thymine in a methanol droplet is summarized in Figure 3. The droplet contained 100 methanol molecules treated at the MM level. The state and orbital characters were the same in the gas phase. Indeed, the S 1 is represented by an n to π * transition while the S 2 state is represented by a π to π * transition. The PES was computed along the geodesic path connecting the ππ * , nπ * and S 0 minimum energy structures as well as the nπ * /ππ * and S 0 /nπ * intersections. The initial bright state in methanol was still the S 2 . From this state, a barrierless access to the nπ * /ππ * intersection was observed, similar to what was observed in the gas phase. From this intersection, a branching towards the nπ * and ππ * minima occurred. From the ππ * state, the S 0 /ππ * intersection was not found. This could be since this intersection was reached by the motion of a CH 3 group (vide supra), which was hindered in solution. Rather, the motions around the ππ * minimum could reach back to the nπ * /ππ * minimum and thus be trapped in the nπ * minimum. From the nπ * minimum, a larger barrier was obtained to reach the S 0 /nπ * than in the case of thymine in the gas phase. The difference between the nπ * minimum and the corresponding intersection with the ground state was around 1.45 eV, thus the reaction was expected to be slower than in the gas phase. To reach the conical intersection, a large reorganization of the solvent molecules was needed in the path connecting these two minimum energy structures. As an estimate, the difference between the purely MM interaction energies among the solvent molecules was 0.1 eV higher in S 0 /nπ * intersection with respect to the nπ * minimum. Still, if the solvent was relaxed along this coordinate, the barrier disappeared, but the S 0 and the nπ * states intersection disappeared. The two states had an energetic gap of around 0.4 eV with the relaxed solvent. This implies that the solvent had a direct participation in the definition of the conical intersection.  The participation of the solvent on the intersections is clearly shown in Figure 4, in which the intersection geometries, branching plane vectors and PES around the intersection points are represented. In methanol, both the S 0 /nπ * and nπ * /ππ * intersections peaked. On the one hand, the nπ * /ππ * intersection was completely planar, unlike the boat-like structure in the gas phase. The participation of the solvent in this intersection was negligible. On the other hand, the S 0 /nπ * intersection was similar to the gas phase, in which the C 4 -O 4a bond was perpendicular to the thymine bond. While the gradient vector of this intersection did not involve the solvent molecules, the approximate non-adiabatic coupling was completely delocalized over all the solvents. This implies that solvent coordinates were to be considered in the definition of the branching plane of conical intersections. Similar conclusions were extracted from model Hamiltonians of conical intersections of protonated Schiff bases in solution by Burghardt et al [44].

Photochemistry of Thymine in Water Droplet
The photochemistry from the lowest energy singlet states of thymine in a water droplet is summarized in Figure 5. The droplet contained 300 water molecules treated at the MM level. The state and orbital characters were the same in the gas phase and methanol droplet. The PES along the geodesic path connecting the ππ * , nπ * and S 0 minimum energy structures as well as the nπ * /ππ * and S 0 /nπ * intersections is shown in the figure. From the bright S 2 state, a barrierless access to the nπ * /ππ * intersection was observed, similar to gas phase and methanol. At variance with methanol and similar to the gas phase, the minimum of the ππ * state was not found in water, indicating that the minimum was close to the S 1 /S 2 intersection. Still, the S 0 /ππ * intersection was not found either, probably for the same reason as in methanol, due to the hindrance of the CH 3 group motion due to the steric interactions with the solvent molecules. From the nπ * /ππ * intersection, all trajectories should go towards the direction of nπ * minimum. A double barrier occurs in the path to the S 0 /nπ * intersection, indicating a larger reorganization of the solvent molecules than in methanol. Indeed, from the nπ * /ππ * to the nπ * minimum a first barrier was obtained, which disappeared when the solvent molecules were relaxed along the geodesic path. Additionally, from the nπ * minimum a larger barrier was obtained to reach the S 0 /nπ * than in the case of thymine in the gas phase. The difference between the nπ * minimum and the corresponding intersection with the ground state was around 1.34 eV, thus the reaction was expected to be slower than in the gas phase, but slightly faster than in methanol. The participation of the solvent on the intersections is clearly shown in Figure 6, in which the intersection geometries, branching plane vectors and PES around the intersection points are represented. In water the intersections were similar to the methanol ones, both the S 0 /nπ * and nπ * /ππ * intersections peaked. On the one hand, the nπ * /ππ * intersection was completely planar, unlike the boat-like structure in the gas phase. The participation of the solvent in this intersection was negligible. On the other hand, the S 0 /nπ * intersection was similar to the gas phase, in which the C 4 -O 4a bond was perpendicular to the thymine bond. While the gradient vector of this intersection did not involve the solvent molecules, the approximate non-adiabatic coupling was completely delocalized over all the solvent, as it was observed from methanol. This implies again that solvent coordinates were to be considered in the definition of the branching plane of conical intersections.

Conclusions
A new development of TD-DFT/MM energy and gradient using electrostatic potential fitting interaction Hamiltonian has been developed and implemented. This method has been used to describe the photochemistry of thymine in the gas phase and embedded in methanol and water nanomeric droplets, using critical point search in the excited states and the potential energy surface cut along a geodesic pathway connecting all minima. We show that the gas phase thymine can undergo an ultrafast relaxation pathway from the S 2 state to the ground state via a series of conical intersections. In solution, the barrierless path to the S 2 /S 1 intersections is similar to in the gas phase. On the contrary, the S 0 /S 1 intersections are highly dependent on the solvent state. Indeed, the approximate non-adiabatic coupling vector defining the branching plane of the S 0 /nπ * intersection strongly mixes the motions of thymine and the solvent. The S 0 /ππ * intersection was not found in solvent, but rather, a transfer from the ππ * minimum to the nπ * minimum is observed.
In conclusion, TD-DFT/MM is an efficient method to describe the excited state potential energy surface of chromophores embedded in complex media. From this study, we inferred that the nπ * /ππ * intersection is independent of the solvent state and therefore should be an equivalently fast reaction both in gas phase and in solution. On the contrary, to reach the nπ * minimum and S 0 /nπ * intersection requires a relaxation of the solvent and a particular configuration of the solvent molecules to reach the intersection. This implies that the S 0 /nπ * intersection will be much slower in solution than in the gas phase. It remains to be seen what implication this has for the photochemistry of nucleobases in DNA. Indeed, the ground/excited intersections involve out-of-plane motions that are strongly hindered in solution, but they might still be possible in the electrostatic environment of DNA. In the future, we will use similar techniques to study the degree of delocalization of the non-adiabatic coupling vectors in DNA.

Data Availability Statement:
All data is available upon reasonable request directly to the author.