To apply our workflow and obtain reference results for benchmark purposes, we carried out calculations on all three solute-solvent systems at the PBE0 56-58 / pcseg-2 53 level of theory. We chose this functional due to its accuracy in the modeling of molecular geometries, and the basis set because it has been shown to yield good results with DFT for both molecular structures and vibrational properties. The grid used was a standard type grid in LSDalton called GRID5, which has a radial integration accuracy of 2.15443 • 10 -17 and an angular expansion order of 47. The inner region had a radius of 16 Å (this choice is based on the results presented in section 4), and the outer region had a radius of 30 Å for all three solute-solvent systems. For the inner region, atom-centered multipoles up to quadrupoles and atom-centered dipole-dipole polarizabilities were derived using the LoProp scheme, based on properties calculated with the Dalton program package employing PBE0 and a recontracted version of 6-31+G* basis set named loprop-6-31+G* in Dalton. The generation of the embedding potentials was facilitated using PyFraME. The Lennard-Jones 12-6 parameters for the inner and outer regions, as well as the charges used only for the outer region, were taken from the GAFF2 force field or TIP3P in the case of water.