Modeling solubility of CO2 gas in room temperature ionic liquids using the COSMOSAC-LANL model: a first principles study

In this paper we present a thermodynamic model for asymmetric solutions with a special emphasis on solute–solvent interactions. The new ‘‘COSMOSAC-LANL’’ activity coefficient model is rooted in first principles calculations based on the COSMO model where the microscopic information passes to the macroscopic world via a dielectric continuum solvation model followed by a post statistical thermodynamic treatment of self-consistent properties of the solute particle. To model the activity coefficient at infinite dilution for the binary mixtures, a 3-suffix Margules (3sM) function is introduced to model asymmetric interactions and, for the combinatorial term, the Staverman–Guggenheim (SG) form is used. The new ‘‘COSMOSAC-LANL’’ activity coefficient model has been used to calculate the solubility of CO2 in room temperature ionic liquids and to model the selectivity between CO2 and CH4 gases. We have shown improved solubility and selectivity prediction of CO2 and CH4 gas in room temperature ionic liquids using the ADF-COSMOSAC-2013 model with the new ‘‘LANL’’ activity coefficient model. The calculated values have been compared with experimental results where they are available.


Introduction
Asymmetric interactions between solute-solvent particles in a solution are responsible for the non-ideal behavior of the solution. 1,2 The chemical potential of species A can be written as m A = m A * + RT ln x A for an ideal solution, whereas for a regular solution (deviation from ideal behavior) the chemical potential is expressed as m A = m A * + RT ln a A . In both equations m A * is the standard state chemical potential and a A is the activity of species A, which is proportional to the mole fraction via the activity coefficient, g A , a A = x A g A . The activity coefficient of a species (g A ) is a measurement of the non-ideality of the solution. Therefore an accurate activity coefficient model is an essential tool to model equilibrium thermodynamical properties, such as solubility, partition coefficients, osmotic pressure, phase equilibrium, and pK a . 3,4 The most popular activity coefficient models in the literature are UNIQUAC (UNIversal QUAsiChemical), 5 UNIFAC (UNIQUAC Functional-group Activity Coefficients), 6 NRTL (non-random two-liquid model), 7 Equation of State (EOS) model, [8][9][10][11][12][13][14][15] COSMO-RS (COnductor like Screening MOdel for Real Solvents) 16,17 and COSMO-SAC (COnductor like Screening MOdel-Segment Activity Coefficient). 18 The first three models are based on the Group Contribution Method (GCM), [5][6][7] while the two COSMO models are based on electronic structure calculations followed by a statistical mechanical treatment of the self-consistent properties of the solute and solvent species. Among these five activity coefficient models, the GCM models are very efficient and fast to compute when the B1000 parameters needed for a chemical compound are known. 19 Based on quantum mechanical calculations, the COSMO-based models are able to separate the isomers and are thus applicable to a large number of systems with fewer parameters. The down side of using the GCM and COSMO models is that they do not predict the strong asymmetric interactions responsible for regular solution behavior found in solutions like CO 2 gas in ionic liquids. In the current version of these models, no explicit term for asymmetric interaction is present and the asymmetric interactions are implicit in the activity coefficient term.
In this paper we introduce a new model that includes a 3-suffix Margules (3sM) 1,20,21 function to explicitly describe asymmetric interactions along with the conventional combinatorial energy term proposed by Staverman-Guggeinheim (SG). 22,23 The former takes into account the interaction energy arising due to the hydrogen bond interactions, dispersion and coulombic interactions between solute and solvent molecules, whereas the entropic effects of size and shape differences present between solute and solvent molecules are included in standard form via the latter term. The Margules parameters were calculated with the COSMOSAC-2013 24 and thus, in combination with the new activity coefficient model, it is denoted as the ''COSMOSAC-LANL'' model throughout this paper. Finally, the proposed model has been applied to calculate the CO 2 solubility in room temperature ionic liquids for which the experimental results were available in the literature. The calculated solubilities using the new model have been compared with other theoretical results in this field.
This article has been organized as follows: the details of theory and computation are reported in Sections 2 and 3, respectively, followed by a discussion on the conformational dependence of the s profile of ionic liquids in Section 4. The results and discussion on the solubility isotherm and s profile are given in Section 5.1 and a discussion on asymmetric behavior in Section 5.2. The solubility results of CO 2 in RTILs are given in Section 5.3. The selectivity of CO 2 over CH 4 and screening of ionic liquids for best CO 2 solubility are given in Sections 5.4 and 5.5, respectively. A brief conclusion is given in Section 6.

Solubility in asymmetric solutions
The solubility of greenhouse gases was recently reported using a new model here referred as the COSMOSAC-LANL 25 activity coefficient model. This model includes two terms: a 3-suffix Margules (3sM) 20,21 function as a quantitative measurement of inherent asymmetric interaction present in a solution due to the hydrogen bond interaction, strong dispersion and coulombic interaction and the long-range Staverman-Guggenheim 22,23 (SG) combinatorial interaction due to the size and shape differences between solute and solvent molecules. We will work in the regime of low solubility where g LANL where we define ln(g asym i/S ) = (DG asym i/S À DG asym i/i )/RT which is the difference between the asymmetric interactions in a mixture (i/S) and in pure state (i/i), which represents the solvation free energy change in terms of solute and solvent interactions when a solute particle goes into a fixed position in solution from a fixed position in its ideal state. Since the ''LANL'' activity coefficient model has asymmetric interaction (ln(g asym i/S )), the total activity coefficient (ln(g LANL i/S )) model is also called an asymmetric model. For pure species, the asymmetric interaction is zero, DG asym i/i = 0, and hence, ln(g asym i/S ) = (DG asym i/S )/RT = (DG solv i/S À DG solv i/i )/RT = (DG solv )/RT. In COSMOSPACE 26,27 (the pictorial representation is shown in Fig. 1), it can be seen that the activity coefficient due to the asymmetric interaction is equal to the activity coefficient due to the short range residual interaction (g residual i/S ) at infinite dilution (x solute=i -0, x solvent=j -1). In COSMOSPACE, the activity coefficient can be written in terms of the segment interaction as where i stands for pure compound, n stands for the different type of segment, and n i stands for the number of segments of type n and say, (ln g n i/S À ln g n i/i ) = lng n i/S (asym). The asymmetric interaction mentioned above can be written in terms of asymmetric segment interaction where DG n i/i (asym)/RT = 0 for an ideal case, so ln g asym i=S ¼ P n n n i DG n ðasymÞ=RT ¼ P n n n i ln g n i=S ðasymÞ À ln g n i=i ðasymÞ .
Since for similar types of segments the concept of asymmetric interaction between segments does not exist, so ln g n i/i (asym) will vanish and therefore the equation will be reduced to ln g asym i=S ¼ P n n n i ln g n i=S ðasymÞ. In this model, we used the 3-suffix Margules function to capture the asymmetric interaction in terms of Margules parameters. Therefore, the expression of the excess Gibbs energy due to the asymmetric interaction can be written as where x 1 and x 2 are the mole fractions for solute and solvent molecules in a binary mixture and A 21 and A 12 are the Margules parameters. For asymmetric solutions when A 12 a A 21 and when they are equal (A 12 = A 21 ), the solution is symmetric reducing eqn (4) to the 2-suffix Margules function. A pictorial representation showing how these Margules parameters are obtained is given in Fig. 1. By differentiating eqn (4) with respect to the mole number for species 1 and 2, 27,28 one can get the asymmetric activity coefficient where n is the total mole number and n i is the mole number of species (i = 1, 2), and x i = mole fraction = n i /n of i-th species. Therefore the expressions for the activity coefficients for both species 1 and 2 present in the binary mixture solution are where a 1 = (2A 21 À A 12 ), b 1 = (2A 12 À 2A 21 ), a 2 = (2A 12 À A 21 ), and b 2 = (2A 21 À 2A 12 ). These parameters are calculated from the activity coefficients of species 1 in species 2 and species 2 in species 1 under infinite dilute conditions. 20 The differentiation of eqn (4) and the derivation of the 2-suffix Margules function from the 3-suffix Margules function are given in Sections 1 and 2, respectively, in the ESI. † Within the 3-suffix Margules law, at infinite dilution eqn (6) and (7) can be written as ln g asym i/S and which is either A 12 (accounts solubility of species 1 in 2) or A 21 (accounts solubility of species 2 in 1). A description has been given in Appendix A.
Hence, one can implicitly express the residual interaction within the asymmetric 3-suffix Margules function for a binary mixture. We use the activity coefficient at infinite dilution (g N ) calculated using the COSMOSAC-2013 24,29 model for all binary mixture solutions used in this work. The reasons for using the COSMOSAC-2013 model are as follows (i) this is the only COSMOSAC model whose parameters are optimized over a large data set for ADF calculation, 24 (ii) this is the only COSMO-SAC model which includes a sophisticated dispersion interaction along with a modified residual and combinatorial term and also, (iii) this model holds an exact expression for the segment chemical potential that satisfies over all boundary conditions, and the resulting equation for the activity coefficient obeys the Gibbs-Duhem relationship to maintain thermodynamic consistency. Therefore, the basic background of the proposed model is based on the COSMOSAC-2013 model in this study. The concept of the COSMOSPACE has been used to explain the implicit definition of the residual interaction within the asymmetric term present in the ''LANL'' activity coefficient model within the COSMO-RS framework. The A 12 and A 21 cases have been shown for a spherical solute and solvent molecules in Fig. 1 which is in a hypothetical COSMOSPACE.
It is already known that the 3-suffix Margules function usually represents a less regular and therefore less asymmetric solution. 1 Since solubility is an entropy driven phenomenon the effect of entropy comes from size and shape differences between the solute and solvent species. To consider the effect due to the different sizes and shapes of the solute and solvent species (long range interaction) the Staverman-Guggenheim (SG) 5,[22][23][24]30 combinatorial term has been used here , where x i is the mole fraction of component i; r i and q i are the normalized volume and surface area parameters for species i and z is the coordination number equal to 10. zq i represents the number of nearest-neighbor sites to one of the solute and solvent molecules. In this model for the combinatorial term only the surface area q i has been normalized with q 0 = 79.53. A detailed description of this SG form has been given in Section 3 in the ESI. † After adding the terms for combinatorial and asymmetric interactions described above, we obtain the expression for the activity coefficient for solute (i) and solvent ( j) species in solution as: and ln g LANL Eqn (9) and (10) At, ln g i/S (comb) o 0, ln g i/S (COSMOSAC-LANL) o ln g i/S (COS-MOSAC-2013) since a similar combinatorial term has been used in both the models and based on the derivation in Appendix A, the effect of the combinatorial term thus has been counted here twice as a consequence of the proposed model. Therefore to overcome the effect due to the combinatorial term, we change eqn (9) and (10) based on our proposed assumption on the low solubility region (gLANL for the solute species, for solvent species, respectively. We have defined eqn (9) and (10) as ln g LANL i or j/S (model) and eqn (12) and (13) as ln g LANL i or j/S (compute) for species 1 and 2, respectively. The COSMO volume and surface for the pure species obtained from the COSMOSAC-2013 model have been used to calculate the combinatorial term. Since we use the same combinatorial term in the model for the asymmetric interaction due to the size and shape differences between the solute and solvent species, we changed the asymmetric contribution coming from the residual interaction by changing the standard state of the g asym i/S from 1 to 2.71 and hence by making the asymmetric contribution slightly more asymmetric towards the less solubility region based on our proposed assumption at the very beginning of this model i.e.; g LANL i/S Z 1. Therefore, we focus on obtaining a more accurate activity coefficient at infinite dilution. A relationship between ln g model i/S and ln g compute i/S has been established by ln g i/S (real) = ln g i/S (compute) = 1 + ln g i/S (model), (14) for the solute species in a binary mixture. The details of the derivation has been given in Appendix B. Eqn (14) is valid when the binary system is approaching x 1 -1; x 2 -0 for a species. The limit mentioned for eqn (14) holds the region x 1 -0; x 2 -1 for a species within it according to the proposed model (Appendix B). The other thermodynamic properties like excess Gibbs energy, free energy of mixing, variation of the logarithm of activity coefficient with the solute mole fraction for the whole range of concentrations obtained using eqn (12) and (13)  with varying mole fractions of the solute species. 1 We calculated the total excess Gibbs energy (G ex ) for the binary system using the equation and the free energy of mixing for the binary mixture using where ln g LANL i/S and ln g LANL j/S are obtained from eqn (12) and (13). At infinite dilution, Henry's constant is related to the infinite dilution activity coefficient and the fugacity of pure gas, 1,2 lim x!0 f i ðT; P; xÞ The solubility of gas molecules in ionic liquids can be measured by taking the inverse of Henry's constant at a partial pressure of 1 bar, i.e., x i = 1/H i/S . Henry's constant is directly related to the free energy of solvation (DG N solv ) which is defined as the change in Gibbs energy upon transferring a gas molecule from the pure gaseous phase at some standard pressure p 0 to the infinite solution and they are related to each other as follows The fugacity of the gas molecule was calculated following the relation described in the work of Lin et al., 31 which is The computational details of the parameters A i , B i , C i , D i and E i are given in the work of Lin et al. 31 This thermodynamic property is a function of solute. The activity coefficients at infinite dilution calculated using the asymmetric model and the fugacity calculated using the above relations have been used to calculate the solubility of CO 2 in ionic liquids. We have calculated the activity coefficients at infinite dilution from the new COSMOSAC-LANL model which has been discussed in detail in the previous section.

Computational details
All the COSMO 19 files for CO 2 , CH 4 and ionic liquids have been generated using the Amsterdam Density Functional software version 2016. 24,32 The initial structures of each gas and ionic liquid have been drawn in ADFview. Then they have been relaxed using a simple UFF method already implemented in ADF2016. Using the final geometry resulting from the UFF calculations, the rest of the calculations have been performed. The geometry optimization was performed using the generalized gradient approximation (GGA) with the BP exchange correlation functional. The zero order regular approximation (ZORA) method was used for calculating the relativistic effect. The TZP basis set for small frozen core and Becke integration with the spline Zlm fit for density fitting have been used in ADF for each molecule followed by the post COSMOSAC calculation already implemented in ADF for the calculation of the s profile of each molecule. It should be noted that the parameters used in the COSMOSAC-2013 24 model are not optimized for ionic liquids. So, in our calculations we simply used those parameters which were optimized for the simple organic molecules and their systems in the ADF software for the COSMOSAC-2013 model. 24 The details of quantum COSMO settings have been given in the paper on the COSMOSAC-2013 model of Xiong et al. 24 Each ionic liquid has been treated as a single intact ionic liquid molecule in the quantum COSMO calculation since it is already known from the work of Kirchner et al. 33,34 that the total charge of each species in an ionic liquid is always less than 1 due to the intermolecular charge transfer between the cation and anion. To ensure this fact, we calculated the s profile for 18 ionic liquids by first treating them as (case I) single intact molecules and then as (case II) separate cation and anion moieties and calculated their s profiles. The s profile (p s ) is defined as the probability of finding a segment of the COSMO surface with charge density s. 35 In most cases, significant differences in their s profile for the different types of treatments (case I) and (case II) of ionic liquids were observed as shown in Fig. 2. The results for case I and case II for 24 ionic liquids are shown in Table 1. The solubility calculations have been done using ADF-COSMOSAC-2013. In case I, each ionic liquid has been treated as a single molecule, therefore their solutions have been treated as a binary mixture. The solution has been treated as a ternary mixture for case II. Therefore, necessary calculations have been done to convert the solubility data by dividing g N by 2. This has been discussed in Section 4 in the ESI. † We notice that the solubility results vary significantly for the different treatments of the ionic liquid using the COSMOSAC-2013 model. The calculated average absolute relative deviations (AARD) for case I and case II are 31% and 63%, respectively. For solubility calculations in ionic liquids, eqn (17) was solved for a partial pressure of 1 bar. In the present calculation, the solubility of two greenhouse gases mainly CO 2 and CH 4 has been calculated at different temperatures and at constant low partial pressure (1 bar) of the gases. All Margules parameters (A 12 ) and (A 21 ) were calculated using the ADF-COSMOSAC-2013 model at x i -0, x j -1 and x j -0, x i -1, for species i and j, respectively. After getting those parameters, we calculated the solubility and other equilibrium thermodynamical properties stated in eqn (15)-(17) (G ex , DG mix and solubility isotherm) for all binary mixtures used in this study. The calculated values of asymmetric interactions along with the calculated combinatorial term at infinite dilution have been given in Table 2. A minimum criterion has been established for when to apply eqn (9) in place of eqn (12) using the relationship in eqn (14) and that is ln g asym i/S (compute) = 1 (Appendix C). Only the ionic liquid ([Eohmim][BF 4 ]) showing the lowest solubility of CO 2 is found to satisfy the above condition.
If a solution obeys Raoult's law, it is known as an ideal solution. Any positive and negative deviation from Raoult's law (p vs. x plot) represents a regular solution. If the deviation from Raoult's law is more, the solution will show asymmetric behavior either in the positive or in the negative direction in the p vs. x plot. Based on the amount of the deviation, it will be known as the symmetric (less deviation from Raoult's law) or asymmetric (more deviation from Raoult's law) regular solution. In the present model, the residual interaction presents implicitly within the asymmetric interaction of the 3-suffix Margules function type. Later, by combining the 3-suffix Margules function with the SG term, we consider all types of short range and long range interactions present in a binary mixture. Such explicit treatment of the asymmetric interaction only between the solute and solvent species allows us to explain any positive or negative deviation (asymmetric interaction) from Raoult's law using the other derivations given in Appendix A to Appendix C in the main article. The 3-suffix Margules function often represents the less regular solution. The less regular behavior of the ADF-COSMOSAC-2013 model has been shown in Table 2. If the ADF-COSMOSAC-2013 model could predict the strong regular solution behavior of the CO 2 /IL binary mixture, then g N (COSMOSAC-2013) a g N (3-suffix Margules function). Since g N (COSMOSAC-2013) = g N (3-suffix Margules function) as shown in Table 2, the COSMOSAC-2013 model presents a less regular solution and thus a less asymmetric solution model for CO 2 solubility in ionic liquids. For thermodynamic consistency it is important to point out that the differentiation must satisfy the Gibbs-Duhem equation The numerical differentiation has been done on a set of activity coefficient data for varying mole fractions for both solute and solvent molecules. An analytical solution of 3-suffix Margules functions and a demonstration of the numerical differentiation for two different systems have been given in Sections 5 and 6, respectively, in the ESI. † We used eqn (12) to calculate the solubility of CH 4 to check the selectivity between CO 2 and CH 4 gases because the asymmetric interactions will be more for CH 4 gas in ILs than the asymmetric interactions observed for CO 2 in ILs in the regime g LANL Z 0. This is indicated from their very high values of fugacity calculated using eqn (19) and lower molecular weights with respect to CO 2 (i.e.; m CO 2 4 m CH 4 ) given in Table 3. The fugacity is inversely related to the solubility of the gas according to eqn (17).

r profile of ionic liquids and their conformational dependence
The impact of different conformations of ionic liquids on the other physical and thermodynamical properties has been looked at and confirmed for 18 different ionic liquids for two different conformers. To check any kind of such dependency, we performed two types of quantum COSMO calculations: (i) the quantum COSMO calculations have been directly performed on  the geometry drawn in ADFview, followed by a relaxation using the UFF method followed by the ADF-COSMOSAC-2013 calculation without any a priori separate geometry optimization and analytical frequency calculation and (ii) the quantum COSMO calculations have been directly performed on the geometry drawn in ADFview, followed by a relaxation using the UFF method followed by the ADF-COSMOSAC-2013 calculation with a priori separate geometry optimization and analytical frequency calculation. The conformation generated following the first method (i) is known as Conf1 and the conformation generated following the second method (ii) is known as Conf2.
Using these schemes, we first did our quantum calculations for geometry optimization followed by the calculations of the analytical frequency to check any presence of the imaginary frequency in the calculated analytical frequency and thus to ensure that the ground state global geometry has been reached. A similar type of quantum settings has been chosen for the calculations because the geometry optimization in the gas phase is already included in the COSMOSAC-2013 model in ADF. Therefore to maintain the parity between these two calculations, the same quantum setting has been selected for the geometry optimization followed by the analytical frequency calculations. The quantum setting for the geometry optimization and analytical frequency calculation have been given in Table 4. The difference between the s profile of the conformations Conf1 and Conf2 (Dp s ) along with the s profile of the conformation Conf2 (p s ) has been shown in Fig. 3. The calculated analytical frequency for 18 different ionic liquids has been shown in Fig. S1 in the ESI. † All 18 ionic liquids studied along with their COSMO surface points have been given in Fig. S2-S4 (ESI †) and the different conformers Conf1 and Conf2 have been shown in Fig. S5 in ESI. † Almost for all cases, we find that the difference between the s profiles (Dp s ) passes through zero without showing any significant changes except for the ionic liquid containing [BF 4 ] À and [PF 6 ] À anion. The calculated solubility and activity coefficient for 18 different ionic liquids for two different conformers have been given in Table 5. From our calculations, it is clear that the conformational dependency on the s profile and hence on the activity coefficient at infinite dilution is very negligible except for the ionic liquid ([Eohmim][BF 4 ]). The reason behind this observation has been explained in Section 5.1. This observation is also reflected on their solubility data given in Table 5 for two different conformers. Since for 18 different ionic liquids, we did not observe any significant changes in the solubility for two different conformers, therefore for the remaining 15 ionic liquids we did our calculations following Conf1 case.
A similar approach to the conformational dependency on the s profile of ionic liquid has been already proposed to model the solubility of metal complexes in ionic liquids using the COSMOSAC-LANL model. 25 In case of the metal complex solubility in ionic liquids, DG ex due to the combinatorial term is significant and higher with respect to the combinatorial interaction between the CO 2 and the ionic liquid. The COSMO volume and surface of a transition metal complex is almost B10 times higher than the same for CO 2 and CH 4 gases. These differences in the COSMO volume and surfaces will effect the combinatorial interaction in those solid-liquid equilibria significantly. We notice that the combinatorial interaction between the CO 2 gas and the ionic liquid has a stabilizing effect while the same between the metal complex and the ionic liquid has a destabilizing effect on the solubility. The result is shown in Table 6.

Solubility isotherm and r profile
The solubility of two gases in room temperature ionic liquids has been calculated using at infinite dilute condition. The all quantities present in the above equation have been explained in eqn (17) in Section 2.
The solubility isotherm has been calculated in all cases by using p i = x i g N i f i (T,P), 31 where g N i is the activity coefficient at infinite dilution, f i is the fugacity and p i is the pressure. The results have been discussed in the following paragraphs.
The solubility isotherms for CO 2 in [C 4 mim][NTf 2 ] ionic liquids at four different temperatures are shown in Fig. 4 along with the experimental results. 36 From the solubility isotherm it is confirmed that the model works up to B6 bar partial pressure very well for CO 2 solubility in the ionic liquid. Similar calculations have been done for the solubility of CH 4 gas molecules in [C 4 mim][PF 6 ] at room temperature (298 K) as shown in Fig. 5. The calculated results are in good agreement with the experimental results at a low partial pressure (1 bar) and hence, in the low solubility region for both CO 2 and CH 4 gas.
In this study we have taken the s profile for each molecule calculated using the COSMO setting mentioned in the work of Xiong et al. 24 Fig. 6 where s is the screening charge density and s 0 is equal to 0.007 e Å À2 for the Gaussian distribution. The relation between them is p s (HB) = p s (OHÁ Á ÁOH) + p s (OHÁ Á ÁOT) and thus p s (total) = p s (HB) + p s (NHB), where p s (NHB) is the sigma profile due to the nonhydrogen bonded group. p s ðHBÞ ¼ A HB i ðsÞ A i P HB ðsÞ and The details are given in the work of Xiong et al. 24 The molecules having the ability to form intermolecular hydrogen bonds between the hydroxyl (OH) group of the cation and the FB group of BF 4 À are the least CO 2 absorbers because the physical absorption of CO 2 molecules in those ionic liquids will cost high Gibbs free energy of cavity formation at the expense of breaking of intermolecular (HÁ Á ÁF) hydrogen bonds. Such kind of (OHÁ Á ÁOT) bond formation has been shown in Fig. 6; plot (a) indicated by a blue line, the smaller amount of (OHÁ Á ÁOH) is indicated by a green line and the total of (OHÁ Á ÁOT) and (OHÁ Á ÁOH) is indicated by a red line and the total s profile p total (s)  showing medium and maximum CO 2 gas solubility have no hydroxyl groups (OH). Therefore all those hydrogen bond interactions in those ILs occur due to the other than hydroxyl group (OT) interactions that are weak and these results are already seen in their corresponding s profiles in Fig. 6; plot (b-d). The qualitative relation between the CO 2 gas solubility and the inverse of the hydrogen bond interaction in an ionic liquid is in good agreement with the other results already present in this field. 39-41

Asymmetric behavior of the solution
In Fig. 7, we have shown the excess Gibbs energy with varying solute mole fractions for three binary mixtures of CO 2 and ionic liquid mixtures showing (1) maximum, (2) medium and (3) least CO 2 gas solubility, respectively. In all three cases, we observed that the excess Gibbs energy shows a minimum at mixture compositions other than 50 : 50 composition of the binary solution and also the free energy curve does not symmetrically decay at the two ends. This asymmetric behavior of the excess free energy vs. x 1 plot is in agreement with the asymmetric behavior of the model. To confirm this asymmetric behavior of the solution, we also plot the logarithm of the   Fig. 7 arises due to the asymmetric nature of the solution, otherwise it would be a straight line. 1 We also plot the logarithm of activity coefficient with x 1 for both the solute and solvent species. The two curves are found to intersect each other at a composition different from the 50 : 50 mixture composition in Fig. 8. All these results for asymmetric interactions indicate towards the different types of interaction present between the solute and solvent molecules in those vapor-liquid equilibria involving CO 2 in ionic liquids.
The new model was used to observe the change in G ex and DG mix as a function of solute concentration. The ionic liquid showing least CO 2 solubility is found to show positive G ex and DG mix over some solute concentration with respect to the ionic liquid showing maximum CO 2 solubility.

Solubility of CO 2 in RTILs at low partial pressure
In this work, the CO 2 solubility in ionic liquid at a low partial pressure (1 bar) of the gas has been calculated. We use our model to calculate the solubility of CO 2 in 33 ionic liquids for which experimental data are available for comparison. The experimental results of solubility have been taken from the references stated in the work of Lin et al. 31 All the calculations have been done in two phases: (i) the solubility at room temperature (293.15 K to 303.15 K) (Fig. 9) and (ii) the solubility for a wide range of temperatures (283 K to 333.15 K) and at an atmospheric pressure of 1 bar (Fig. 10). Fig. 9 and 10 show a comparison of COSMOSAC-LANL with experimental data and the previous model COSMOSAC-2013. The % errors in the linear regression for COSMOSAC-LANL and COSMOSAC-2013 are 17.2% and 30%, respectively. The percentage error is calculated using the relation Fig. 10, it is also clear that the predictive and accuracy power of the COSMOSAC-LANL model is good to produce the experimental solubility with respect to the COSMOSAC-2013 model. Table 6 shows the results obtained using the COSMOSAC-LANL and COSMOSAC-2013 models for the solubility of CO 2 in ionic liquids at room temperature and 1 bar partial pressure. In the case of COSMOSAC-LANL and COSMOSAC-2013 models, calculations on solubility and Henry's constant at room temperature and over a wide range of temperatures have been performed. Accuracy of the model has been tested by calculating the error percentage (average absolute relative deviation (AARD)),    (17), we know that x i = 1/H i/S (T,P) = 1/(g N i/S f i (T,P)). Therefore, Henry's constant is directly related to the activity coefficient at infinite dilution at a particular temperature. Fig. 11 shows the direct correlation between experimental solubilities and 1/g N i/S for both COSMOSAC-LANL and COSMOSAC-2013 models at a particular temperature and at 1 bar partial pressure. A good correlation has been observed using the new COSMOSAC-LANL model (R 2 = 0.88). We also calculate the % error in the linear regression for the ideal solvation case when g CO 2 = 1 for the two models at an average temperature of 298. 15  ] ionic liquids, for the rest of the ionic liquids improved solubility results were obtained using the COSMOSAC-LANL activity coefficient model. The possible reason for the deviation in those cases is the combinatorial interaction present in the model. For those four ionic liquids, the results obtained using the COSMOSAC-2013 model are in good agreement with the experimental results. We found that the model can be modified for those four ionic liquids by considering only 3-suffix Margules expression without any SG term present in the model for the asymmetric interaction according to the data in Table 2. We observed that the calculated solubility at room temperature and at 1 bar partial pressure is found to improve to 11% for 29 from 13% for 33 selected ionic liquids. The AARD for solubility is Fig. 8 (a-d)   improved to 13% for 29 from 14% for 33 selected ionic liquids for a wide range of temperatures and at 1 bar partial pressure.

Ideal-dilute selectivity of CO 2 over CH 4
The study of absorption selectivity between the CO 2 over the CH 4 gas is important to large scale industries such as the natural gas sweetening process involving the selectivity between CO 2 and CH 4 . 31, 42 We used our model for calculating the selectivity of CO 2 /CH 4 for which the experimental data are available for comparison at a low partial pressure of 1 bar. The selectivity has been calculated from the single solubility data of a gas and thus by taking the ratio of that solubility data, Fig. 10 Experimental vs. calculated solubility of CO 2 in different ionic liquids at varying temperatures (283 K to 333.15 K) and at 1 bar partial pressure.   Table 9. All the experimental results on selectivity have been taken from the references stated in the work of Lin et al. 31 The present method is found to predict the selectivity between CO 2 and CH 4 within 12% error at a low partial pressure of 1 bar and in the low solubility region.

Screening of ionic liquids based on the solvophobic property
We proposed a new method to screen ionic liquids for the best solubility of CO 2 gas capture based on the solvophobic effect of the ionic liquids reported in the work of Sedov et al. 43 Solvophobic effect is the generalized concept of hydrophobic effect for self associating ionic liquids and the key driving force of self-assembly of amphiphilic compounds into the mesophase structures such as micelles, vesicles, bilayers, microemulsions and emulsions. 44 Similar to the recent work of Sedov et al. 43 in 2016, we define the solvophobic effect in this work. In ionic liquids, the concentration of ions (the cation plus the anion of the liquids), C i = 1/V m , which is the inverse of molar volume of the ionic liquid, is correlated with thermodynamic solvation properties (DG solv ). Since the Henry constant and free energy of solvation are related to each other following the thermodynamic relation given in eqn (18), the experimental free energy of solvation of CO 2 for all ionic liquids has been calculated and correlated to the inverse of the molar volume of the ionic liquids in Fig. 12. We have also shown the correlation between the experimental solubility and the inverse of the molar volume Fig . 11 Experimental solubility of CO 2 in different ionic liquids as a function of (1/g N i/S ) at 298.15 K at 1 bar partial pressure.  Fig. 13. A correlation between the calculated and the experimental free energy of solvation has been shown in Fig. S6 in the ESI. † The calculated free energy of solvation has been obtained using the both COSMOSAC model in this study. The % error in the linear regression is 0.43% in the COSMOSAC-LANL model and 21.45% in COSMOSAC-2013 (the result is shown Fig. S6, ESI †) for free energy solvation. The molar volume has been calculated using the following equation V m (cm 3 mol À1 ) = V COSMO (Å 3 per molecule) Â (10 À24 cm 3 Å À3 ) Â 6.02 Â 10 23 molecules per mol. The calculated COSMO volume, temperature, experimental solubility and experimental free energy of solvation have been given in Table 10 separately. We compared the experimental DG solv (exp.) vs. 1/V m (exp.) and DG solv (exp.) vs. 1/V m (cal.) for 16 different room temperature ionic liquids for which the experimental densities were available under ambient conditions. We also plotted 1/V m (exp.) vs. 1/V m (cal.) and calculated the AARD. The results are given in Fig. S7(a, b) and Table S3 in the ESI. † We found that the COSMOSAC model underestimates the 1/V m (cal.) by 14.1%. The plot shown in Fig. S7(b) in the ESI, † indicates small improvement in the results of the inverse of experimental and calculated molar volume. The model shows significant improvement when we compared the experimental and the calculated DG solv obtained from both the COSMOSAC-2013 and COSMOSAC-LANL models. The result is shown in Fig. S6 (ESI †). The COSMOSAC-LANL model predicts the free energy of solvation within 4% error while the COSMOSAC-2013 model predicts it within 7% error. Also the correlation between the experimental and the calculated free energy of solvation improves for the COSMOSAC-LANL model (R 2 = 0.87) almost by 10% with respect to the COSMOSAC-2013 model (R 2 = 0.78). The molar volume has been calculated using the scheme mentioned in the link of COSMO-RS tutorial available in SCM (https://www.scm.com/doc/ Tutorials/COSMO-RS) for ionic liquids, 32 however we made a little change in our calculation for single molecule treatment of the ionic liquid. The molecular weight and density of 16 different ionic liquids has been taken from the IL Thermo database. 45 We have screened the ionic liquids for CO 2 gas for the ionic liquids mainly containing [NTf 2 ] À anion, [C 2 mim] + and [C 4 mim] + cations for which the experimental solubility data were available. We observed that these are the few common cations and anions usually used for CO 2 gas capture using ionic liquids. A relationship  has been established between the inverse of the calculated molar volume (1/V m ) (x) obtained from the quantum COSMO calculations and the experimental solubility data (y) for each case. A correlation has been observed between these two properties for all the ions and results are given in ESI, † in Fig. S8-S10. The linear relationships to screen the cations (mainly imidazolium and pyrrolinidium) against a [NTf 2 ] À anion is y = À0.007x + 0.0571 with R 2 = 0.77. The same to screen the different anions for [C 4 mim] + and [C 2 mim] + cations are y = À0.0062x + 0.0541 for R 2 = 0.96 and y = À0.0057x + 0.052 for R 2 = 0.90, respectively. Therefore, one will be able to compute the solubility of an unknown ionic liquid having these cation or anion common in them using these relations, provided the molar volume is calculated using quantum COSMO calculation within the ADF-COSMOSAC-2013 model implemented in ADF-2016 software. Also one can predict the solubility of CO 2 gas in any ionic liquid very well according to the relation (y = À0.0062x + 0.0536 with R 2 = 0.82) shown in Fig. 13. We observe that the outlier present in the Fig. 13 is the [C 6 tma][NTf 2 ] ionic liquid which is a dication dianion ionic liquid. The coefficient of determination is improved by 10% (R 2 = 0.92; result is not shown here) when we assume the molar volume/ cation-anion pair present in such type of ionic liquid.

Conclusion
A new thermodynamic model is presented to express the solution properties for CO 2 and CH 4 in room temperature ionic liquids (RTILs) under varying thermodynamic conditions. The new thermodynamic model is composed of the Staverman-Guggenheim (SG) term due to the long-range interaction and short range asymmetric interaction (3-suffix Margules function) between the solute and solvent species in a solution. For the short range asymmetric interaction, we chose the 3-suffix Margules function for the asymmetric solution because we assume that the long range combinatorial term will be responsible for the entropy driven solubility phenomenon due to the size and shape differences between the solute and solvent species. We have shown in this study that the new thermodynamic model for activity coefficient improves the solubility prediction of CO 2 in ionic liquids when it was used with the ADF-COSMOSAC-2013 model. We calculated the solubility from the Henry constant available from the fugacity and activity coefficient at infinite dilution for the two greenhouse gases mainly CO 2 and CH 4 in RTILs. We take 33 such example cases to test our model which is comprised of a long range combinatorial term along with the short range asymmetric interaction. Our present model was not only able to improve the quantitative prediction of solubility of CO 2 in those 33 ionic liquids, but was also able to explain other thermodynamical and physical properties of those systems. It has been found from our study that the present model improves the solubility of CO 2 gas by 17% with respect to the COSMOSAC-2013 model at room temperature over 33 points. For the quantitive prediction of Henry's constant, the calculated percentage of errors is (1) 14% for COSMOSAC-LANL and (2) 22% for COSMOSAC-2013 at room temperature over 33 points, respectively. For a wide range of temperatures, the present model predicts solubility within 14% AARD with respect to the COSMOSAC-2013 model over 79 points used in this study. With this new model the selectivity between CO 2 /CH 4 was predicted within 12% error over 18 points. We have proposed a new scheme to screen ionic liquids for the best CO 2 solubility using the concept of solvophobic effect which is a correlation between the solubility and the inverse of the molar volume of the ionic liquid resulting from the first principles calculation in this model in addition to the conventional screening method which is a correlation between the solubility and the 1/g N i/S at a particular temperature. Using the new COSMOSAC-LANL thermodynamic model the solvation mechanism of CO 2 and CH 4 in those ionic liquids has been explained. The calculated theoretical results are found to be in good agreement with the other experimental and theoretical results where they were available. The only drawback of this model is that it cannot be applied to model the systems tending towards the supercritical condition because the present model within the COSMO-RS frame work cannot handle the free volume effect and therefore the vapor-liquid transition. 46 However, the model can be applied to a certain range of pressures and temperatures for both CO 2 and CH 4 gases to predict their solubility in ionic liquids.
The new thermodynamic model can be applied to study the other vapor-liquid, liquid-liquid and solid-liquid equilibria. Very recently, the new COSMOSAC-LANL model has been applied to predict the solubility of metal complexes in ionic liquids to screen the metal complexes 25 and ionic liquids for the best solubility of the redox active species in ionic liquids. The solvation mechanism of those systems has been explained using the new COSMOSAC-LANL model.

Conflicts of interest
There are no conflicts to declare. and ln g model 1 = ln g model 1 (comb) + ln g model 1 (asym).
Since we did not change the term due to combinatorial contribution, ln g compute 1 (comb) = ln g 1 model (comb).