advanced the field of computational chemistry, they do exhibit certain limitations that stem from three important facts. First, they primarily consist of molecules that are considerably smaller than what is commonly encountered in modern medicinal chemistry. Second, their structures have been optimized using theoretical models that do not account for molecule-solvent and collective dispersion interactions. Lastly, they have not fully explored the vast conformational landscape inherent in these molecules. Especially, the interaction between the molecule and the chemical environment (i.e., solvent) is crucial when investigating molecules of pharmaceutically relevant size, as it is well-known that drug binding does not occur in vacuo but rather in solutio. Omitting solvent effects can result in inappropriate treatment of conformations, tautomers, or molecular reactivity. In the computational modeling of molecular systems, solvents can either be considered implicitly or represented explicitly. However, owing to the computational cost and intricate nature of explicitly representing the solvent, most QM studies opt for the utilization of implicit solvent models such as conductor-like screening model for real solvents (COSMO-RS) , modified Poisson-Boltzmann (MPB) , and Generalized Born (GB) model augmented with the hydrophobic solvent accessible surface area term (GBSA) .