Examining the corresponding reorganization energies in the bottom part of Table , we see that the highfrequency reorganization energies, λ HF , are around 0.05 eV for all MM solvent environments, suggesting that the contribution from the higher-frequency chromophore vibrational degrees of freedom is consistent across all solvent environments. The low-frequency reorganization energy, λ LF , from the CV + chromophore only without any MM solvent, is only 0.006 eV, reflecting a minimal contribution from chromophore vibrations to the lowfrequency spectral density. Including the four peripheral H-bonding solvent molecules barely changes the reorganization energies, indicating that the strong peripheral hydrogen bonds contribute little to the low-frequency solvent modes for inhomogenous broadening, consistent with the simulated optical spectra. On the other hand, the inclusion of the four axial solvent molecules increases λ LF to 0.030 eV, close to the value for full MM solvent (0.040 eV). The λ LF value for the single H-bonding axial solvent molecule is 0.024 eV, close to that from four axial solvent molecules and not far from the total MM solvent λ LF value, showing that a single axial solvent molecule is responsible for well over half of the λ LF reorganization energy of a full solvated system. This result once again shows that the low-frequency solvent modes, specifically the motions of the axial solvent molecules that presumably interact strongly with the π-electrons of the chromophore, cause the long-time relaxation of the energy-gap fluctuation and the resulting inhomogeneous broadening in the absorption spectrum. It is well-known that in the slow modulation (inhomogeneous) limit, the Stokes shift, as defined by the difference between absorption and emission maxima, is twice the reorganization energy. Therefore, we estimate that the solvent-induced Stokes shift from our QM/MM AIMD simulation would be about 2λ LF = 0.080 eV, not far from the observed experimental Stokes shift of 0.084 eV 89 or 0.099 eV for CV + in methanol. Overall, our results demonstrate that it is critical to accurately account for the interactions between the chromophore and the axial solvent molecules in both MD simulations and excited state calculations in order to properly capture solvent-induced inhomogeneous broadening in the spectrum.