Molecular modeling was performed at the density functional theory (DFT) level using consistently the ORCA 4.2.1 code. Geometry optimization was performed for the S0, S1 and T1 state at density functional theory (DFT) level for S0 and T1, or its Time-Dependent extension in the Tamm-Dancoff Approximation (TDA) for S1. Two microhydrated tautomeric forms of the main chromophore involving a network of two additional water molecules have been considered. The wB97x exchange-correlation functional including Grimme D3 dispersion corrections and the 6-311+G(d,p) basis set has been consistently used. The vibrational frequencies have been obtained within the harmonic approximation by diagonalizing the hessian matrix, calculated at the same level of theory as for the geometry optimization. Non-adiabatic molecular dynamic was performed in vacuo within the Surface Hopping (SH) formalisms using the SHARC code coupled with ORCA. To assure statistically significant sampling and to simulate the effects of a nuclear wave-packet an ensemble of 200 independent non-adiabatic trajectories has been considered. The initial conditions have been obtained sampling the Franck-Condon region of the non-tautomerized M in presence of two water molecules through a Wigner distribution. 50 non-adiabatic dynamic trajectories have been propagated computing potential energies and gradients at TDA level, explicitly calculating the SOC element at each step, with the wB97x/6-31G level of theory, with a time-step of 0.5 fs for a total time of 2 ps (4000 steps). Persico decoherence corrections have been consistently applied. Note that, as shown in SI, the reduced level of theory has been chosen since it provides a correct energetic ordering of the states and a proper description of the involved PES, as compared to the higher one. For the SH nonadiabatic dynamic propagation, our chosen model explicitly comprised the first excited singlet state, and the lowest two triplet states. An initial full population of the S1 state was enforced, and the time evolution of the population of the different states has been fitted with exponential functions to obtain the corresponding kinetic model and the characteristic life-times.