## 2 Materials and Methods  
### COSMO Simulation  
A total of 1298 solvent combinations composed of type 5 HBA (total of 22) and HBD (total of 59) at a 1:1 ratio were used in this study (Table S1). The HBAs consist of quaternary ammonium salts, quaternary phosphine salts, terpen, and imidazolium salts. The HBDs consist of alcohols, organic acids, and profens. These components all exhibit significant van der Waals-capable surfaces which give them the potential to form a hydrophobic DES complex. The sequence of computations was performed using the COSMOLogic platform in series: TurbomoleX (TmoleX19), COSMOConfx (COSMOConf18), and COSMOThermX. The molecules that were not available in the COSMOtherm onboard database were input to TmoleX19 in the SMILES format, as found from PubChem. The lowest geometrical energies were computed along with sigma surfaces at the bp86 functional and Karlsruhe (Ahlrichs) def2-TZVP (default-2 Valence Triple-Zeta Polarization) basis set, as recommended in literature [28]. The output files (.cosmo extension) were input to COSMOConfx for conformer analysis at the same theory levels. The results were then uploaded to COSMOThermX for property predictions of the studied systems.  
Equation (1) was used to compute chemical potentials \((\mu_{\rm s}(\sigma))\) of the DES-difluoromethane complex by means of sigma surface interactions, which would be used as a basis for successive property predictions.  
\[\mu_{\rm s}(\sigma)=-\frac{\rm RT}{a_{\rm eff}}\ln[\int P_{\rm s}\big{(}\;\sigma^ {\prime}\big{)}\mathrm{e}^{\big{(}\frac{\mathrm{d}\sigma}{\rm RT}\;(\mu_{\rm s }\;(\;\sigma^{\prime})-E_{\rm mish}(\sigma,\;\sigma^{\prime})-E_{\rm HB}(\sigma,\;\sigma^{\prime}))\big{)}}\mathrm{d}\sigma^{\prime}] \tag{1}\]  
where \(E_{\rm mish}\) is the misfit energy, hydrogen bond energy is \(E_{\rm HB}\), and van der Waals energy is \(E_{\rm vdw}\); \(\sigma\) and \(\sigma^{\prime}\) are two interacting surface segments between two molecules x and x\({}^{\prime}\), and \(p_{\rm s}(\sigma)\) is the distribution function or sigma profile. R is the ideal gas constant and T is absorption temperature. \(a_{\rm eff}\) is the area of effectiveness term.  
The activity coefficient \(\big{(}\gamma_{\rm S}^{\rm i}\big{)}\) was calculated through the difference in the chemical potentials of component difluoromethane in solvent (s) and pure form (p), as seen in Equation (2).  
\[\ln\gamma_{\rm s}^{\rm R32}=\frac{\Big{(}\mu_{\rm s}^{\rm R32}-\mu_{\rm P}^{ \rm R32}\Big{)}}{\rm RT} \tag{2}\]  
In Equation (3), [PERSON]'s constant (H) was used for the validation procedure and was calculated through an iterative process of varying the pressures (P) with the concentration (C).  
\[\mathrm{C}=\mathrm{H}\times\mathrm{P} \tag{3}\]  
The predicted solubility of R-32 in DES is produced through Equation (4).  
\[x_{\rm R32}=\frac{P_{\rm R32}}{P_{\rm R32}^{\rm o}\times\;\gamma_{\rm R32}} \tag{4}\]  
where \(p_{\rm j}\) is the partial pressure of R-32 in the DES system, \(p_{\rm j}^{\rm o}\) is the partial pressure of R-32 in its pure form, \(x_{\rm j}\) is the solubility of R-32 in HDES, and \(\gamma_{\rm j}\) is the activity coefficient of R-32 in solvent.  
Finally, the excess enthalpy of mixing H\({}_{\rm mix}\) was determined through the addition of three enthalpic contributors available in Equation (5).  
\[H_{\rm mix}=H_{\rm mf}+H_{\rm hb}+H_{\rm vdw} \tag{5}\]  
here, these contributions are from enthalpic change due to the misfit of sigma charge segments H\({}_{\rm mf}\), enthalpic change due to hydrogen bonding (H\({}_{\rm hb}\)), and contribution due to van der Waals interactions (H\({}_{\rm vdw}\)).