Multiscale methods, such as combined QM/MM (quantum mechanical/molecular mechanical) methods, are used to model a part of the entire system at the QM level of theory while the rest is treated using a classical MM approach. In this way, QM/MM methods provide a robust description of localized quantum phenomena in complex environments with feasible computational effort. A lot of effort currently goes into the development of multiscale QM/MM methods both to extend their robustness and applicability. In the presented work, we resort to the QM/MM model with the polarizable embedding, which accounts for the mutual polarization between the QM and MM subsystems. This type of description is essential to account for the charge density response of the QM subsystem due to the electrostatic changes in the environment. Such effects play a critical role in processes like absorption and emission and can affect the properties of molecules with charge separation, etc. The development of polarizable embedding models 4 continued over the last few decades. Initial developments based on polarizable continuum models were then modified and adapted for explicit, heterogeneous environments represented by Drude oscillators, fluctuating charge, or point induced-dipole models. Although such a sophisticated representation of the environment in hybrid QM/MM methods comes with a higher computational cost, the continuous increase in available computational power also fuels interest in more precise models.