A Plausible Mechanism of Uracil Photohydration Involves an Unusual Intermediate

It is well-known that photolysis of pyrimidine nucleobases, such as uracil, in an aqueous environment results in the formation of hydrate as one of the main products. Although several hypotheses regarding photohydration have been proposed in the past, e.g., the zwitterionic and “hot” ground-state mechanisms, its detailed mechanism remains elusive. Here, theoretical nonadiabatic simulations of the uracil photodynamics reveal the formation of a highly energetic but kinetically stable intermediate that features a half-chair puckered pyrimidine ring and a strongly twisted intracyclic double bond. The existence and the kinetic stability of the intermediate are confirmed by a variety of computational chemistry methods. According to the simulations, the unusual intermediate is mainly formed almost immediately (∼50–200 fs) upon photoabsorption and survives long enough to engage in a hydration reaction with a neighboring water. A plausible mechanism of uracil photohydration is proposed on the basis of the modeling of nucleophilic insertion of water into the twisted double bond of the intermediate.

All energies are relative to the respective ground state minimum energy and are given in electron volt (eV). See Figure

Remarks about performance of MRSF-TDDFT
The mixed-reference spin-flip (MRSF) TDDFT enables accurate computation of the excitation energies S9 and the real (and avoided) crossings between the potential energy surfaces of the ground and excited states, the so-called conical intersections, S10 of organic molecules. Simultaneously, MRSF-TDDFT eliminates the major drawback of the spin-flip TDDFT methods, the spin-contamination of the ground and excited states, without loosing the efficiency and transparency of the formalism. S11,S12 MRSF has the ability to accurately reproduce the relative energy (with respect to the Franck-Condon energy) of conical intersections, as shown in comparison between the MRSF and MRCISD+Q energies in Ref. S10 Especially valuable for studying excited state dynamics is the ability of MRSF to correctly describe the adiabatic potential energy surfaces in the vicinity of the Franck-Condon geometry and far from it. S13-S15 In a series of recent studies, S16,S17 the accuracy of MRSF was benchmarked against the results from three independent types of the equation-of-motion coupled-cluster S-10

Analysis of the optimized structures
The S 0 equilibrium geometry of uracil obtained with MRSF-TDDFT/BH&HLYP/6-31G* method features planar structure, which is in excellent agreement with the experimental crystallographic geometry and high level theoretical calculation; see Figure S2 and Table. S1. As the consequence of populating the π * 1 orbital, the C 4 =C 5 bond is partially broken and this leads to lengthening of the C 4 =C 5 bond (See Figure S1 for frontier orbitals of uracil and Figure S2 for optimized structures and atom numbering). The neighboring C 5 -C 6 and C 4 =O 8 bonds are simultaneously shortened and lengthened as the π * 1 orbital has bonding and anti-bonding characteristic between these atoms respectively. This unique skeletal vibration is denoted as bond length alternation (BLA) and defined as where ∆R's are displacements with respect to the S 0 equilibrium geometry. Along the BLA motion, conical intersection (CI) between the S 2 and S 1 states (CI 21, BLA ) and excited state S 1 minimum can be found. They have flat ring with a BLA distortion of 0.13 and 0.21 respectively. Therefore, CI 21, BLA is close to Franck-Condon (FC) geometry while S 1 minimum (S 1, min ) is relatively far from FC geometry and exhibits an inversion of double bond(C 4 =O 8 , C 5 =C 6 ) and single bond(C 4 -C 5 ). It has been recently shown by some of us S17 that CI 21, BLA is essentially involved in the very early dynamics of thymine. Figure S1: Frontier orbitals of Uracil (n O7 , π 1 , π 2 , and π * 1 ).

S-11
Meanwhile, hydrogen out of plane (HOOP) torsion about the C5=C6 bond denoted as φ coordinates produces a CI 21, ± with an angle of ±48°. As both are the mirror images of one another, in the following, only the structures with positive values, e.g. CI 21, + , will be considered. The CI between S 1 and S 0 state (CI 10, + ) and ground state local minimum (S 0, min + ) are also found with a more HOOP distortion of φ = 118.21°and 176.3°respectively.
Recent work by Chakraborty et al. S1 who recently employed extended multi-state CASPT2 (XMS-CASPT2) in the non-adiabatic molecular dynamics (NAMD) simulation showed that the ethylenic CI 21 has a very small deviation from planarity, while CI 10 is considerably non-plnar.  The φ values are given in parenthesis next to C 5 =C 6 bond in (°).

Vertical excitation energies
Extensive experimental and theoretical works are collected in Table S2. The three lowest excited states (S 1 , S 2 , and S 3 ) are characterized by n → π * , π 1 → π * 1 , and π 2 → π * 1 transitions respectively (See Figure S1 for frontier orbitals of uracil). The failure of CASSCF in the description of uracil photo-dynamics is originated from the overestimation of S 2 state due to the lack of dynamic electron correlation, resulting in S 2 -trapping model. S1 This is because the bright S 2 (π 1 → π * 1 ) state has a strong ionic character, requiring a proper description of dynamic electron correlation. However, CASSCF is especially deficient for this, which inevitably overestimate the S 2 -S 1 energy gap (∆ 12 ). As expected, the CASSCF produce overestimated values of 1.63 ∼ 1.64 eV (See Table S2) as compared to the experimental value of 0.70 eV, S18 S18 it should be remembered that the experimental estimate of the n → π * transition corresponds to the origin of the transition band, but the theoretical values correspond to vertical transitions. Therefore, the vertical gap ∆ 12 from the experiment may be somewhat smaller than 0.7 eV; thus, improving the agreement between theory and experiment. Taking into account our previous NAMD study of thymine and its consistency with the high level ab initio computations, MRSF-TDDFT in connection with the BH&HLYP density functional was selected for the NAMD simulations in the present work. S17 S-16 Table S2: Vertical excitation energies of the three lowest singlet (n O7 π * 1 , π 1 π * 1 , π 2 π * 1 ) excited states of uracil obtained with various methods from the literature. All energies are relative to the respective ground state minimum energy and are given in electron volt (eV). See Figure S1 for frontier orbitals of uracil.
S-17 Table S3: Relative energies of the three lowest singlet (S 2 , S 1 , S 0 ) states of uracil obtained by MRSF with different basis set and functional. The relative energies (in eV) are given relative to the respective ground state energy in the BH&HLYP/6-31G* optimized geometries.   Figure S5: Minimum energy paths (MEPs) on the S 2 (blue), S 1 (red), and S 0 (black) PESs of uracil obtained using the nudged elastic band (NEB) optimization S2,S3 in connection with the MRSF-TDDFT/BH&HLYP/6-31G* method.

Results of MRSF-TDDFT gas phase NAMD Simulations
The four-state (S 0 , S 1 , S 2 , and S 3 ) NAMD simulations were initiated by sampling the Wigner function of a canonical ensemble S26,S27 at T = 300K around the S 0 equilibrium geometry. The nonadiabatic coupling vectors are computed numerically by using a fast overlap calculation. S28 The NAMD has been performed based on the fewest-switches surface-hopping algorithm. S29  S-20 CASPT2 (XMS-CASPT2) starts with ca. 40 % of S 2 state, the MRSF-TDDFT starts with ca.
90 % of S 2 state, which perhaps makes the τ 1 (the S 2 → S 1 internal conversion process) of the former (12.5 fs) slightly smaller than the latter (37.0 ±1.2 fs) (See Table. S4). In the MRSF-TDDFT NAMD simulations, within the first ca. 20 fs, the population of the S 2 state drops to around 40 % (due to the S 2 → S 1 population transfer) and becomes almost negligible at t 100 fs. In both simulations, the S 0 state population starts to appear at ca.
40 fs and MRSF-TDDFT remarkably well reproduce the overall dynamics; the rise of S 0 , decrease of S 1 and even their population cross between 150 ∼ 250 fs of XMS-CASPT2 as shown in Figure S6.
Regarding the ground state population changes, as discussed in the main text,

Generality of occurrence of 4
Although the formation of an intermediate similar to 4 has not previously been observed in NAMD simulations of other pyrimidine bases, e.g., thymine and cytosine, S17,S25,S44 its absence among the products obtained theoretically does not necessarily imply the impossibility of its existence. Indeed, typical duration of the NAMD simulations is   Figure S10: Geometries of torsionally distorted structures 4 (4-ura) of uracil, 4-thy of thymine, and 4-cyt of cytosine. The key bondlength (inÅ), the diherdal angle θ (see Figure 1d of the main article for definition), and the relative energy (∆E in kcal/mol) with respect to the ground state equilibrium structure are given.
Nevertheless, the structures similar to 4 can, in principle, be obtained in the geometry optimizations started with a suitable geometry. If the optimized structure is S-28 a local minimum, such geometry optimization can be a proof of the existence of the intermediate 4 for other pyrimidine nucleobases. To verify this conjecture, a series of geometry optimizations have been performed with the BH&HLYP/cc-pVTZ method for torsionally distorted structures of thymine and cytosine. These optimizations did indeed produce structures similar to 4 for the two other pyrimidine bases. The structures shown in Figure S10 are the local minima on the ground electronic state PES; they were characterized by vibrational analysis, which did not produce imaginary frequencies.
The 4-thy and 4-cyt structures occur at approximately the same relative energy with respect to the ground state equilibrium geometry as 4. Hence, these calculations confirm feasibility of existence of the torsionally strained intermediates in other pyrimidine bases. Figure S11 shows the IR spectra of uracil and intermediate 4. The spectra were built using the experimental matrix isolation IR frequencies and intensities from Ref. S5 and the IR frequencies (cm −1 ) and intensities (km/mol) computed by the BH&HLYP/cc-pVTZ method in the gas phase. The computed frequencies were scaled by the factor 0.9335. S6 The spectra were convoluted with the Lorentzian lineshape function with half-width at half maximum 10 cm −1 .

S-46
Cytosine and Thymine