The functional reactions in blue light photoreceptors generally involve transiently reduced flavins exhibiting characteristic infrared (IR) spectra. To approach a theoretical understanding, here we apply density functional theory (DFT) to flavin radicals embedded in a molecular mechanics (MM) model of an aqueous solution. Combining a DFT/MM approach with instantaneous normal-mode analyses (INMA), we compute the IR solution spectra of anionic and neutral flavin radicals. For a set of mid-IR marker bands, we identify those changes of spectral locations, intensities, and widths, which are caused by sequentially adding an electron and a proton to the oxidized flavin. Comparisons with experimental IR solution spectra of flavin radicals show the accuracy of our DFT/MM-INMA approach and allow us to assign the observed bands. The room temperature ensembles of solvent cages required for the INMA calculations of the IR spectra are generated in an MM setting from molecular dynamics (MD) simulations. For the solvated flavin radicals, these MD simulations employ MM force fields derived from DFT/MM calculations.