We then analyzed the respective Dyson orbitals 42 to assign the computed IEs to particular water molecules. Fig. shows a structure of 5 water molecules embedded in the MM part with the Dyson orbital on the central water molecule. The assignment of IEs can also be accomplished by considering the Mulliken charges from the natural orbital analysis 62 (we used this approach in the production simulations). We then constructed the bulk spectra by taking the IEs corresponding to a specified number of water molecules. For example, in the calculations with 5 QM waters, one can construct the spectrum by taking all 5 IEs or by taking only the IE of the central water molecule. The overall spectra were constructed as histograms (with 0.05 eV bins) by collecting the IEs from the snapshots and then convoluted with gaussians. The spectra shown in the manuscript were produced using gaussians with 0.2 eV width (FWHM). Discarding some of the computed IEs in the construction of the spectra slows down the convergence with respect to the number of snapshots, but removes the artifacts due to waters on the QM/MM interface. As we show below, the best protocol (in terms of balancing accuracy of IEs and convergence with respect to sampling) entails using 20 QM waters and assembling the spectra by taking the IEs of the 5 central waters.