The HF/pcseg-1 39 level of theory was used for the full QM calculations on the small model systems. For the QM/MM calculations on the large systems and the small model systems, the solutes were also modeled at the HF/pcseg-1 level of theory, while the effects from the surrounding solvent environment were included using the fragment-based PE model. Additionally, for formaldehyde in water, we also considered a system where the two water molecules hydrogen-bonded to the formaldehyde are included in the QM subsystem and thus also treated at HF/pcseg-1 level of theory. The embedding potentials representing the solvents were obtained using the PyFraME Python package. Each solvent molecule is described by atom-centered multipoles, up to and including quadrupoles, atom-centered dipole-dipole polarizabilities, and a Lennard-Jones 6-12 potential. The multipoles and polarizabilities were derived using the LoProp scheme based on PBE0/6-31+G* calculations performed with the Dalton program package. The Lennard-Jones parameters were sourced from the GAFF force field, as provided by Van der Spoel et al. on , with the exception of water, for which the CHARMM TIP3P Lennard-Jones parameters were applied. The geomeTRIC package was used together with the LSDalton program for the geometry optimizations. The environmental contributions for the QM/MM geometry optimizations were included using the FraME library. The small model systems were geometry optimized using full QM and QM/MM. Note that the solvent was frozen in all QM/MM optimizations.