The solubility of caffeine in aqueous binary mixtures was measured in five aprotic proton acceptor solvents (APAS) including dimethyl sulfoxide, dimethylformamide, 1,4-dioxane, acetonitrile, and acetone. The whole range of concentrations was studied in four temperatures between 25 °C and 40 °C. All systems exhibit a strong cosolvency effect resulting in non-monotonous solubility trends with changes of the mixture composition and showing the highest solubility at unimolar proportions of organic solvent and water. The observed solubility trends were interpreted based on the values of caffeine affinities toward homo- and hetero-molecular pairs formation, determined on an advanced quantum chemistry level including electron correlation and correction for vibrational zero-point energy. It was found that caffeine can act as a donor in pairs formation with all considered aprotic solvents using the hydrogen atom attached to the carbon in the imidazole ring. The computed values of Gibbs free energies of intermolecular pairs formation were further utilized for exploring the possibility of using them as potential solubility prognostics. A semi-quantitative relationship (R2 = 0.78) between caffeine affinities and the measured solubility values was found, which was used for screening for new greener solvents. Based on the values of the environmental index (EI), four morpholine analogs were considered and corresponding caffeine affinities were computed. It was found that the same solute–solvent structural motif stabilizes hetero-molecular pairs suggesting their potential applicability as greener replacers of traditional aprotic proton acceptor solvents. This hypothesis was confirmed by additional caffeine solubility measurements in 4-formylmorpholine. This solvent happened to be even more efficient compared to DMSO and the obtained solubility profile follows the cosolvency pattern observed for other aprotic proton acceptor solvents.