The photochemistry of nucleobases, important for their role as building blocks of DNA, is largely affected by the electrostatic environment in which they are soaked. For example, despite the numerous studies of thymine in solution and DNA, there is still a debate on the photochemical deactivation pathways after UV absorption. Many theoretical models are oversimplified due to the lack of computationally accurate and efficient electronic structure methodologies that capture excited state electron correlation effects when nucleobases are embedded in large electrostatic media. Here, we combine mixed-reference spin-flip time-dependent density functional theory (MRSF-TDDFT) with electrostatic embedding QM/MM using electrostatic potential fittingfitted (ESPF) atomic charges, as a strategy to accurately and efficiently describe the electronic structure of chromophores polarized by an electrostatic medium. In particular, we develop analytic expressions for the energy and gradient of MRSF/MM based on the ESPF coupling using atom-centered grids and total charge conservation. We apply this methodology to the study of solvation effects on thymine photochemistry in water and thymine dimers in DNA. In the former, the combination of trajectory surface hopping (TSH) nonadiabatic molecular dynamics (NAMD) with MRSF/MM remarkably revealed accelerated deactivation decay pathways, which is consistent with the experimental decay time of ∼400 fs. The enhanced hopping rate can be explained by the preferential stabilization of corresponding conical interactions due to their increased dipole moments. Structurally, it is a consequence of characteristic methyl puckered geometries near the conical intersection region. For the thymine dimer in B-DNA, we found new photochemical pathways through conical intersections that could explain the formation of cyclobutadiene dimers and 6−4 photoproducts.