All COSMO-RS calculations were performed using the COSMOtherm software (version 19.0.4, with the BP_TZVP_19.ctd parameterization, COSMOlogic, Leverkusen, Germany). To begin with, the quantum chemical Gaussian09 package was employed to optimize the structures of the studied compounds, which include CO₂, N₂, and components of IL, at the B3LYP/6-31++G (d, p) level. Frequency calculations were conducted to confirm that the optimized structures correspond to true minima. Second, the resulting COSMO files of the optimized structures were subsequently imported into the COSMOtherm program to compute the solubility of CO₂ and N₂ in the studied ILs. For the solubility calculations of gases in ILs, the cation and anion components are treated as separate molecules with equal molar fractions (n_cation = n_anion = n_IL), furthermore, the input variables (T, P) were set to be consistent with the experimental conditions reported in the literature.

At present, multiple ML algorithms have been used to estimate the physical and thermodynamic properties of ILs and IL-involved systems. Among them, the XGBoost algorithm proposed by Chen and Guestrin (2016) is a powerful and efficient algorithm owing to its high training efficiency, good prediction effect, multi-controllable parameters, and user-friendly features. XGBoost can be regarded as a variant of Gradient Boosting Decision Tree (GBDT). Unlike GBDT, XGBoost introduces regular terms to limit the model complexity to reduce the probability of over-fitting, and the second-order derivative information is used for optimization, which accelerates the convergence process of the model and improves the training efficiency. By assuming a dataset contains n examples and m features, the mathematic expressions (Equation 1) and objective function (Equation 2) of the XGBoost algorithm are outlined as follows: y ^ i = ∑ k = 1 K f k x i , f k ∈ F (1)

Here, f _k is the kth independent tree, and F represents the space of regression trees. o b j = ∑ i = 1 n l y ^ i , y i + ∑ k = 1 N Ω f k (2) where l is a differentiable convex loss function that measures the difference between the prediction y ^ i and the target y i , and Ω is the regularization term.

Since the selection and number of features (i.e., the functional groups) significantly affect the accuracy and generalization ability of the ML model, the division of the functional groups followed the JR method in this work (Nannoolal et al., 2007), with the detailed information provided in Supplementary Table S1. Also, for the studied ILs, the same functional group may be contained in both cations and anions, and to better describe the impact of functional groups in anions and cations on solubility, a “-” sign was added after the functional groups from anions. Consequently, the studied ILs were divided into 41 groups for CO₂ solubility modeling and 38 groups for N₂ solubility modeling.

Before model development, the data used were normalized and standardized to eliminate the effects of data magnitude. First, the CO₂ and N₂ solubility datasets were divided into the training set and the test set, with a division ratio of 8:2. The input features for the XGBoost-GC model include temperature (T), pressure (P), and groups on cations and anions (41 for CO₂ dataset, and 38 for N₂ dataset). The target variable for the CO₂ dataset was set to be the relative deviation ( x Exp − x C O S M O − R S x Exp ) between the experimental results and the predictions generated by the original COSMO-RS model. For N₂, the target variable is the absolute deviation ( x Exp − x C O S M O − R S ) between the experimental values and the COSMO-RS model predictions for each sample. For comparison, a model with the same input features but using experimental values as the target variables was also studied, which is named XGBoost-GC-D.

The optimal parameters were obtained through the Bayesian optimization algorithm. Since the XGBoost algorithm is a decision tree-based model, the number of trees should be proper, and too few trees will result in poor prediction, while too many trees may lead to over-learning and over-fitting. The same goes for the maximum depth of the tree. Therefore, simultaneous optimization was performed on these parameters, where the number of trees ranged from 1 to 100 with corresponding maximum depth from 1 to 10. The ranges for the learning rate and subsample ratio were set to 0.01–0.3 and 0–1, respectively. The number of iterations was 200, and the specific parameters are listed in Supplementary Table S2.

Given that the work of Lei et al. (2014) systematically collected CO₂ solubility data in ILs published before 2013 and used it as a database, the CO₂ solubility data used in this work mainly come from literature reported in the past decade (Nonthanasin et al., 2014; Tagiuri et al., 2014; Gonzalez-Miquel et al., 2014; Makino et al., 2014a; Bahadur et al., 2015; Carvalho et al., 2014; Makino et al., 2014b; Zhou et al., 2014; Zeng et al., 2015; Almantariotis et al., 2017; Liu et al., 2016; Zhou et al., 2016; Zoubeik et al., 2016; Nematpour et al., 2016; Watanabe et al., 2016; Zubeir et al., 2016a; Zubeir et al., 2016b; Dai et al., 2017; Jalili et al., 2017; Bai et al., 2017; Mirzaei et al., 2018; Zhao et al., 2018; Jalili et al., 2019; Nath and Henni, 2020; Safarov et al., 2019; Wang et al., 2022; Henni et al., 2023; Kodama et al., 2023; Mirzaei et al., 2023; Suzuki et al., 2024a; Suzuki et al., 2024b), and the experimental data with zero or negative solubility are not considered. Finally, 3,036 sets of CO₂ solubility (mole fraction: 0.00116–0.713) in 72 different ILs were selected at temperatures of 273.15–413.15 K and pressures of 9.7–6,532.8 kPa.

However, for N₂, the relevant experimental data are much less abundant than for CO₂. Here, we collected and screened N₂ solubility data in the previous literature (Zhou et al., 2014; Almantariotis et al., 2017; Liu et al., 2016; Zhou et al., 2016; Jacquemin et al., 2006a; Jacquemin et al., 2006b; Zhou et al., 2013; Almantariotis et al., 2012; Stevanovic and Gomes, 2013; Zhao et al., 2011; Anderson et al., 2007; Yuan et al., 2006; Afzal et al., 2015; Zhang et al., 2017; Bentley et al., 2023). Similarly, the datasets with zero or negative solubility were discarded. A total of 457 N₂ solubility data points in 31 types of ILs were collected, with values ranging from 0.000171 to 0.6187 mol fraction at 283.20–353.20 K and 4.69–14982 kPa. Supplementary Tables S3, S4 provided the detailed experimental ranges of temperature, pressure, and solubility for various CO₂-IL and N₂-IL systems.

Supplementary Figures S1, S2 show the temperature, pressure, and solubility distributions. It could be seen that the temperature data distribution of the two datasets is relatively uniform, while the pressure data were mainly concentrated in 0–1,000 kPa. The CO₂ solubility data is relatively evenly distributed, while the N₂ solubility data is mainly concentrated below 0.05.

The chemical structures of the cations and anions investigated in this work are illustrated in Supplementary Table S5. The cations include imidazolium, pyridinium, pyrrolidinium, ammonium, and phosphonium, and the anions contain acetate, sulfate, sulfonate, tetrafluoroborate [BF₄], hexafluorophosphate [PF₆], Bis [(trifluoromethyl)sulfonyl]azanide [NTf₂], etc.

Appropriate model evaluation metrics are crucial for evaluating the accuracy of the model. To provide a reasonable evaluation, the average absolute relative deviation (AARD, Equation 3) and coefficient of determination (R2, Equation 4) were used to quantify the discrepancies between the experimental and predicted CO₂ solubilities, where the former is a bias-centric metric while the latter is a variance-oriented one. However, for the N₂ dataset, considering the low accuracy of experimental measurements linked to the low solubility of N₂ in the solvents, the verage absolute deviation (AAD, Equation 5) and R² were used. A A R D % = 1 N ∑ i = 1 N x i − x i ′ x i × 100 (3) R 2 = 1 − ∑ i = 1 N x i ′ − x i 2 ∑ i = 1 N x ¯ i − x i 2 (4) A A D % = 1 N ∑ i = 1 N x i − x i ′ × 100 (5) where N is the total number of samples, the experimental and predicted values of gas solubility in ILs are denoted as x i and x i ′ , respectively, and x ¯ i represents the mean value of the gas solubility in ILs.

As described in Section 2.1, the solubility of CO₂ and N₂ in the identified ILs under the same conditions (T, P, ILs) as reported in the literature was predicted using COSMOthermX (version 19.0.4) and compared with the experimental values (see Supplementary Tables S6, S7). Figures 1A, B present the comparison of the experimentally determined and COSMO-RS predicted gas solubility of CO₂ and N₂, respectively. In Figure 1A, it is evident that the COSMO-RS model tends to underpredict the solubility of CO₂ in ILs, with an AARD of 43.4% and a R² of 0.599. For N₂, as depicted in Figure 1B, the solubility data are spread on either side of the diagonal, with an AAD of 4.95% and a R² of 0.242.

It should be emphasized that the overall trend of the solubilities predicted by COSMO-RS is consistent with the experimental data at various temperatures and pressures (Supplementary Figures S3, S4), confirming that the qualitative prediction of COSMO-RS can be used to reliably screen ILs based on their gas solubilities. To further analyze the model predictions, the N₂ solubility in two ILs was taken as an example to discuss the effects of temperature and pressure. As depicted in Figure 1C, the model prediction performance for [HMIM][eFAP] depends on the studied temperature. As the temperature increase, the points on the consistency diagram get closer to the diagonal line, i.e., the model prediction gets close to the experimental results, and thus the corresponding AAD gradually decreases. This indicates that the prediction of COSMO-RS is more accurate at relatively high temperatures. The same trend was observed for [MDEA][Cl] (Figures 1C, D). Additionally, when the temperature remains constant (e.g., T = 303.4 K), as the pressure increases, both the experimentally measured and theoretically predicted solubilities of N₂ in [HMIM][eFAP] show the same increasing trend and the accuracy of COSMO-RS is gradually decreasing (Supplementary Figure S5). The results of this study demonstrate that, within a certain temperature and pressure range, COSMO-RS can accurately capture the effects of input variables (T, P) on the solubility of N₂ in different cation and anion combinations.

Furthermore, the results of the COSMO-RS model developed in this study were compared with other predictive models reported in the literature. For example, Kamgar and Rahimpour (2016) used UNIQUAC and quantum models to predict the solubility of CO₂ in seven ILs. The study found that UNIQUAC showed good prediction ability for the ILs studied, the ARD in most cases lower than 5%, and the maximum ARD is 9.17%. The predictions of the UNIQUAC model in the literature perform better. Additionally, the COSMO-RS model was also used to predict the CO₂ solubility for the same system, showing an ARD ranging from 6.1% to 62.4%, especially, when the pressure increases, the error becomes larger. Recently, Chen and co-workers used a hierarchical extension strategy to develop a UNIFAC-IL-Gas model for gas solubility prediction. The results showed that for 13 types of gases, including CO₂ and N₂, its prediction performance exceeded the COSMO-RS model (Chen et al., 2020c). The above results further confirm that compared to models that require parameters obtained from the fitting of experimental data, the COSMO-RS model without requirements of any experimental information predicts results qualitatively.

As mentioned before, many studies have demonstrated that higher accuracy can be achieved by performing linear regression on the predicted values obtained by COSMO-RS. These corrected models typically use the experimental values as the target variables. However, in this work, it is evidenced that there is no simple linear relationship between T, P, x C O S M O − R S and x Exp (Supplementary Figures S6, S7), a polynomial expression (Equation 6) combined with different regression strategies were used to further improve the prediction of COSMO-RS: ∆ x = f T , P (6)

For CO₂, the relative deviation ( ∆ x 1 = x Exp − x C O S M O − R S x Exp ) was used (Equation 7): ∆ x 1 = k 1 T + k 2 P + k 3 T 2 + k 4 T P + k 5 P 2 (7)

For N₂, the absolute deviation ( ∆ x 2 = x Exp − x C O S M O − R S ) was used (Equation 8): ∆ x 2 = k 1 T + k 2 P + k 3 T 2 + k 4 T P + k 5 P 2 (8)

Here, k₁-k₅ are the adjustable parameters.

Based on the collected data point, the adjustable parameters were obtained, as listed in Supplementary Table S8. Figure 2 shows the comparison between the experimental gas solubilities and those predicted by the two models. It can be evident from Figure 2A that the CO₂ solubility predictions from the modified model align more closely with the experimental values than those from the original COSMO-RS model. After the model modification, its AARD was decreased to 11.9% with a R² of 0.970. For comparison, the AARD for the original COSMO-RS is 43.4%. For N₂, the modified model shows only a very slight decrease in AAD, and there is no noticeable improvement in R² compared with that before the modification. These results demonstrate that the corrected model improves the accuracy of the COSMO-RS model for predicting CO₂ solubility. However, such a correction does not work for the solubility of N₂ in ILs. The reasons for the above phenomenon are summarized as follows: 1) The CO₂ dataset and the N₂ dataset may have different quality levels. The data in the CO₂ dataset is more accurate and complete and thus can be corrected for better accuracy. 2) The model assumptions themselves and the selection of features are less applicable to the N₂ dataset than to the CO₂ dataset. 3) The insufficient number of samples in the N₂ dataset prevents the model from effectively learning the relationship between the initial predicted value and the experimental value.

The COSMO-RS model can be used for qualitative prediction, which is sufficient for IL screening. The correction with a polynomial expression on COSMO-RS can improve the prediction capability in the solubility for certain gases (CO₂, etc.) but not for all (e.g., N₂). In this section, an alternative option was used to develop a hybrid model, where XGBoost-GC was coupled with COSMO-RS to achieve reliable predictions of CO₂ and N₂ solubility in ILs.

Machine learning has demonstrated significant potential in predicting various properties of ILs, particularly in fields such as green chemistry and electrochemical processes. ILs possess a variety of tunable properties, which are often time-consuming and costly to determine experimentally. ML models, trained on experimental data or theoretical predictions, offer a rapid and efficient means of predicting key properties such as viscosity, density, conductivity, and solubility. However, the performance of ML models is highly dependent on the quality and comprehensiveness of the datasets used for training, and thus the availability of high-quality data remains a critical challenge.

In addition, various thermodynamic models have shown high prediction accuracy for IL-containing systems due to their solid thermodynamic foundations. Effectively combining ML algorithms with these models to improve prediction accuracy without relying on large amounts of experimental data is crucial yet highly challenging.

The accuracy of ML models is highly depended on the selection of meaningful features, such as temperature, pressure, and structural information. The selection of features that better represent the geometric and electronic structures of ILs, along with the application of data-cleaning techniques, can further improve prediction accuracy. Additionally, future advancements may involve the implementation of more sophisticated algorithms, such as deep neural networks, which have the potential to capture complex, non-linear relationships between the structures of ILs and their corresponding properties.