We present a novel framework to self-consistently model the effects of radiation fields, dust physics and molecular chemistry (H$_2$) in the interstellar medium (ISM) of galaxies. The model combines a state-of-the-art radiation hydrodynamics module with a non-equilibrium thermochemistry module that accounts for H$_2$ coupled to a realistic dust formation and destruction model, all integrated into the new stellar feedback framework SMUGGLE. We test this model on high-resolution isolated Milky-Way (MW) simulations. We show that photoheating from young stars makes stellar feedback more efficient, but this effect is quite modest in low gas surface density galaxies like the MW. The multi-phase structure of the ISM, however, is highly dependent on the strength of the interstellar radiation field. We are also able to predict the distribution of H$_2$, that allow us to match the molecular Kennicutt-Schmidt (KS) relation, without calibrating for it. We show that the dust distribution is a complex function of density, temperature and ionization state of the gas which cannot be reproduced by simple scaling relations often used in the literature. Our model is only able to match the observed dust temperature distribution if radiation from the old stellar population is considered, implying that these stars have a non-negligible contribution to dust heating in the ISM. Our state-of-the-art model is well-suited for performing next generation cosmological galaxy formation simulations, which will be able to predict a wide range of resolved ($sim 10$ pc) properties of galaxies.