While molecular dynamics (MD) simulations with explicit solvent molecules (and if needed enhanced sampling techniques) constitute a powerful tool to estimate the conformational ensemble of small molecules in solution , they are inherently constrained by their relatively high computational demand. The molecule of interest (the solute) is solvated in a large box of solvent molecules, and the whole system is propagated through time. This approach is inherently inefficient since the probabilities of the conformational ensemble are identified from their repeated occurrence in the simulation rather than their intrinsic free energy, which is not readily accessible from explicit-solvent simulations. Nevertheless, such simulations are currently the gold standard because implicit solvation methods (e.g, Poisson-Boltzmann (PB) , fast analytical continuum treatments of solvation (FACTS) , or generalized Born (GB) models ) do not achieve the desired accuracy.