In recent years, a lot of effort has been devoted to development and benchmarking of computational procedures aimed at calculating molecular response properties of chromophores in their natural environments such as in a solvent, a protein or another bio-structure. Because molecular response properties are coupled to the electronic structure of the molecule in question, such calculations require at least a partial quantum mechanical treatment, which in many cases is challenging due to the generally large size of the considered molecular systems and the associated computational cost of such calculations.