In summary, through our benchmarked workflow, we provide guidelines for how to use a focused PE-QM approach to simulate harmonic IR and Raman spectra of solute-solvent systems. The benchmark calculations quantify the errors introduced through various boundary conditions imposed by our scheme, such as pseudo-translational and pseudo-rotational contributions, basis-set size, system size, and configurational sample size. Lastly, we compared our theoretical spectra to experimental data, revealing more in-depth strengths and weaknesses. Some aspects warrant further attention. One is the size of the core region when partitioning the system. In this regard, it is unclear what the impact is from the removal of temperature effects and by how much the accuracy can be improved. Of particular interest is how including different types of hydrogen bonds between the core region and the environment region will affect the calculated results.