Solubility model of metal complex in ionic liquids from first principle calculations

A predictive model based on first principles calculations has been proposed to study the solid–liquid equilibria comprising of metal complexes and ionic liquids. The model is based on first principle COSMO calculation followed by post statistical thermodynamical treatment of self-consistent properties of solute and solvent molecules. The metal complex and ionic liquid have been treated as a simple binary mixture. The ionic liquid has been treated here as a single intact molecule. The experimentally observed dual-solute relationship between the ionic liquid and redox active species in presence of a third organic solvent has been established using our model in this work. Within the model, the dual-solute relationship appeared as a simple Gibbs–Duhem relationship between these two species at ambient condition. The dual-solute relationship between the metal complex (V(acac)3, Cr(acac)3 and Mn(acac)3) and ionic liquid ([Tea][BF4]) has been validated by calculating the Gibbs–Duhem relationship, xsolutevs. xsolvent(IL) and 1/γsolutevs. xsolvent(IL) plots. The present model has been applied to a set of ionic liquids, metal complexes and organic solvent (acetonitrile) for which experimental study has been done. The solvation mechanism of the metal complexes in those ionic liquids was obtained using the model. According to our findings, the ionic liquid containing imidazolium cation and [NTf2]− anion is appeared as a suitable solvent for the non-aqueous redox flow cell. We have compared our results with the already reported experimental results where they were available for the non-aqueous solvents.


Introduction
Ionic liquids have emerged as potential non-aqueous media in redox ow batteries.  A key property affecting the performance of a ow cell is the solubility of the redox active material in the supporting electrolytes, therefore, a good solubility model is essential to predict the performance of these systems. Solubility has been modeled in several ways using Hansen solubility parameters, Scatchard-Hildebrand solubility, 22 and rst principle-based COSMO-RS 23-25 calculations followed by statistical mechanical treatment of the self-consistent properties obtained from the quantum calculations for solid dissolution in liquid. However the rst two empirical models depend on parameters to be extracted from experimental data and they depend on the nature of the solvating media. Therefore, if experimental results are not available for a specic solvating media, the solubility of that material cannot be predicted. Also, to predict the solubility of solid in liquids, one needs to have apriori knowledge of heat of fusion and melting temperature of those material, which has limited the use of COSMO 24,25 based solubility models. In this work we show how one can predict a good trend of solubility of the redox active materials in ionic liquids without having any apriori knowledge of heat of fusion and melting temperature using the rst-principle COSMOSAC model. 26,27 These calculations are based on rst-principle approach via a dielectric continuum model and thus have fewer adjustable parameters than the previous models that can be easily adaptable to the metal complex and ionic liquids. These calculations are fast and with good accuracy. Also the solubility measurement of metal complexes in ionic liquids is a challenging task to do experimentally due to the slow diffusivity of the redox active species through ionic liquids in certain cases. 12,16,18 Therefore, the model can be used as a very efficient predictive tool for solubility calculations of metal complex in ionic liquids, to screen ionic liquids for the best solubility of certain metal complexes and to study the solvation mechanism of the metal complexes in the ionic liquid and organic solvent.
In a non-aqueous redox ow cell, the redox active species and supporting electrolyte such as ionic liquids are dissolved in the non-aqueous organic solvent medium. The solubility of the redox active species is affected by the presence of the second solute in the medium because its presence changes the activity of the solution. The presence of this second solute may also affect the mobilities of charged species and hence, ionic conductivity of the solution. According to previous studies, 1,9 the solubility of redox active species decreases with increasing concentration of the supporting electrolyte (IL) which serves as the second solute and the acetonitrile (ACN) as solvent. This phenomenon is known as the dual-solute effect in solubility. 1,9 However, establishing such relationship experimentally is very difficult. Such relationship can be studied directly by measuring the interaction between the redox active species and the ionic liquids in terms of the non-ideality correction (activity coefficient) at innite dilution (g N ). In this paper, we repeat calculations of vanadium tris acetylacetonato (V(acac) 3 ), chromium tris acetylacetonato (Cr(acac) 3 ) and manganese tris acetylacetonato (Mn(acac) 3 ) in non-aqueous medium such as ionic liquids and organic solvent.
In Section 2, the details of solubility theory has been discussed. The computational details have been discussed in Section 3. The discussion on metal complex s prole and on ionic liquids s prole have been given in Subsection 4.1 and 4.2, respectively. The dependence of s prole on the metal complex conformation has been given in Subsection 4.3. The calculations of solubility, screening of ionic liquid and tests of sensitivity with respect to the model parameters have been discussed in Section 5. The theoretical modeling of dual-solute effect has been explained in Section 6. A brief conclusion is given in Section 7.

Theory
The solubility of the metal complexes was calculated using a new model here referred as COSMOSAC-LANL activity coefficient model. This model includes two terms: A 3-suffix Margules (3sM) 28,29 function as a quantitative measurement of inherent asymmetric interaction present in a solution due to the hydrogen bond interactions, strong dispersion, short-range entropic contribution and coulombic interactions and the long-range Staverman-Guggenheim 30,31 (SG) combinatorial interactions due to the size and shape differences between a solute and solvent molecules. We will work in the regime where g LANL i/S $ 1 i.e.; ln g LANL i/S $ 0 since we are dealing with the solid-liquid equilibria at low pressure. Hence, we invoke a new asymmetric interaction in terms of "LANL" activity coefficient model where, we dene ln(g asym i/S ) ¼ (DG asym i/S À DG asym i/i )/RT which is the difference between the asymmetric interactions in mixture (i/S) and 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 xed position in solution from a xed position in its ideal state. Since "LANL" activity coef-cient model has the asymmetric interaction (ln(g asym i/S )), the total activity coefficient (ln(g LANL i/S )) model will be 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. In COSMOSPACE, 32 the activity coefficient due to the asymmetric interaction is equal to the activity coefficient due to the short range residual interaction (g residual ) at innite 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 the pure compound, n stands for the segment n i , and n i stands for the number of segments and say, (ln g i/S n À ln g i/i n ) ¼ ln g i/S n (asym). The asymmetric interaction mentioned above can be written in terms of asymmetric segment interaction where, DG i/i n (asym)/RT ¼ 0 for ideal case, so ln g asym i=S ¼ P n n i n DG n ðasymÞ=RT ¼ P n n i n ðln g i=S n ðasymÞ À ln g i=i n ðasymÞÞ: Since for similar type of segments the concept of asymmetric interaction between segments does not exist, so it will vanish and therefore the equation will be reduced to ln g asym i=S ¼ P n n i n ln g i=S n ð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 free 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. The solution is asymmetric when A 12 s A 21 and when they are equal (A 12 ¼ A 21 ), the solution is symmetric reducing to the 2-suffix Margules function. A pictorial representation about obtaining these Margules parameters has been given in Fig. 1. By differentiating eqn (4) with respect to the mole number for species 1 and 2, 33,34 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 ith species. One can then get the expression for activity coefficients for both species 1 and 2 present in the binary mixture solution and they are where a 1 ¼ ( 12 ). These parameters are obtained from the activity coefficients of species 1 in species 2 and species 2 in species 1 at innite dilute condition. 28 Within the 3-suffix Margules law, at innite dilution eqn (6) and (7) can be written as ln g asym i/S either A 12 (accounts solubility of species 1 in 2) or A 21 (accounts solubility of species 2 in 1), therefore ln g Hence, in COS-MOSPACE, one can show that the asymmetric interaction can be well explained using the residual interaction of COSMO-RS 35 or COSMOSAC 26 and they are equal to the activity coefficient at innite dilution (g N ) and can be called within the 3-suffix Margules function for a binary mixture. We use the activity coefficient at innite dilution (g N ) calculated using COSMOSAC-2013 (ref. 26) model for all binary mixture solutions used in this work. The same Staverman-Guggenheim (SG) 30,31 combinatorial term used in the work of Xiong et al. 26,36 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. The zq i represents the number of nearest-neighbor sites to one of the solute and solvent molecules. In this model for combinatorial term only the surface area q i has been normalized with q 0 ¼ 79.53 because r 0 will be canceled out in the above equation. Aer adding the energy terms for combinatorial and asymmetric interactions described above, we obtain the expression for the activity coef-cient for species i in solution S and they are for solute (i) and solvent (j) species; and ln g LANL The eqn (9) and (10) are called as "LANL" activity coefficient model in the rest of the paper and in some cases the model is referred as COSMOSAC-LANL model as we imported Margules parameters from COSMOSAC-2013 model. At innite dilution ((x solute¼i / 0, x solvent¼j / 1) the two models COSMOSAC-LANL and COSMOSAC-2013 are connected by ln g(COSMOSAC-LANL) À ln g(COSMOSAC-2013) ¼ ln g(comb).
The details of the derivation has been given in Appendix A. It is noteworthy in this regard that in the COSMOSAC-LANL model the asymmetric interaction has been dened explicitly in terms of 3sM function and the residual interaction has been called implicitly within the 3sM function. To verify the asymmetric solution model, we computed various thermodynamical properties (such as G ex and DG mix ) and the correlations between ln g i , ln g j and ln g i g j ! with varying mole fraction of solute species. 22 We calculated the total excess Gibbs free energy (G ex ) of the binary system from the equation (12) and the free energy of mixing for the binary mixture using where x i and x j are the mole fractions of solute and solvent molecules in a binary mixture and ln g COSMOSAC-LANL i/S and ln g COSMOSAC-LANL j/S are the activity coefficients of solute and solvent molecules, respectively obtained from the eqn (9) and (10). The solubility has been calculated using the following equation 37 where, x mc is the mole fraction of the metal complex in the liquid state, DH f mc is the heat of fusion and T m is the melting temperature of the metal complex. For a particular metal complex DH mc and T m are xed values, therefore, if we vary solvent, the solubility is a function of activity coefficient of the metal complex in liquid,

Computational details
All the COSMO 25 les for the metal complexes and ionic liquids were generated using the Amsterdam Density functional soware version 2016. 26,38 The geometry optimization was performed using generalized gradient approximation (GGA) with BP exchange correlation functional. The zero order regular approximation (ZORA) method was used for calculating the relativistic effect. TZP basis set for small frozen core and Becke integration with spline Zlm t for density tting have been used in ADF for each molecule followed by the post COSMOSAC calculation already implemented in ADF for the calculation of sigma prole of each molecule. It is to be noted that the parameters used in the COSMOSAC-2013 (ref. 26) model are not optimized for transition metals and hence for their complexes and for ionic liquids. So, in our calculations we used the parameters which were optimized for the simple organic molecules and their systems in the ADF soware for COSMOSAC-2013 model. The details of quantum COSMO settings have been given in the paper on COSMOSAC-2013 model of Xiong et al. 26 The unrestricted calculation considering spin polarization has been done only for the transition metal complexes as a new addition in the already stated COSMO setting in the work of Xiong et al. All ionic liquids have been treated as a single intact ionic liquid molecule during quantum calculation since the total charge of each cation and anion is always less than 1 due to the intermolecular charge transfer between cation and anion reported in the work of Kirchner et al. 39,40 According to their studies, ionic liquid with different clusters sizes of unity charges is erroneous because signicant charge transfer between the ionic liquid ions reduces their total charge. To double check, we calculated the s prole for 10 ionic liquids by rst treating them as (case I) single intact molecule and then as (case II) separate cation and anion moieties and calculated their s proles. s prole (p s ) is dened as the probability distribution for nding a segment of the COSMO surface with charge density s. 41 In most cases, signicant differences in their s prole were observed for the different type of treatments (case I) and (case II) of ionic liquids shown in Fig. 8 in Section 4.2. Notice that the rst step in the COSMOSAC calculation is geometry optimization. For each COSMOSAC calculation, to ensure the minimization in optimized geometry, we also calculated the analytical frequency to check any presence of imaginary frequency in the calculated spectra for the 11 ionic liquids used in this study. The optimized geometry aer this calculation has been used for the COSMOSAC calculation implemented in ADF for the version 2013. 26 The calculated analytical frequencies for the 10 room temperature ionic liquids for which we propose the solubility model have been given in  Table 1.
The initial conformers of V(acac) 3 ( Fig. 10) metal complex were selected randomly, except the Conf 5 (taken from CSD 42 ), followed by pre-optimization using UFF method and geometry optimization has been done for each conformer followed by the COSMOSAC calculation.

s prole of metal complexes
In spite of the lack of COSMOSAC parameters for transition metals in the metal complex molecules of this study one can safely assume that the ligand will play the major role in the solubility of metal complex because it is the ligand the one interacting with the solvent (see Fig. 2 and 3). The effect of the metal is assumed to be via its effect on the electronic structure of the ligands. Therefore, we used a group contribution method [43][44][45] that assumes that the total s prole of the metal complex is the summation of the s prole of each ligand present in the metal complex. For our study we follow two families of metal complexes: (a) with same ligands but different metal centers and (b) with different ligands with same metal center. In the rst case we studied vanadium, chromium and manganese with acetylacetonate (acac) ligand. In the second case the ligands were derivatives of acac with systematic Ph group substitution coordinating chromium. The structure of the metal complexes and the ligand induced COSMO charge density on the metal center for six different metal complexes studied in this work have been shown in Fig. 3. We also plot the relative differences between the s proles of the metal complexes with same type of ligands. In spite of the similar type of s proles, we observe a signicant difference in their Dp s vs. s plot in Fig. 4. The calculated COSMO surface, volume and molar mass for all metal complexes from I to VI have been given in Table 2. They can be directly correlated with the solubility of a metal complex.
The reason for higher solubility in a metal complex can be explained using the s prole which also present the total COSMO surface of the metal complex. We observed that the metal complexes of almost similar COSMO volumes and surfaces have similar type of s proles while the metal complexes with different type of ligands have different s proles and their total COSMO volume and surfaces are found to increase with the increase of Ph ligand substitution. As seen in Fig. 5, the metal complex with more bulky ligands (VI) is found to show the highest peak height in its s prole which also corresponds to the increasing total COSMO surface. In COS-MOSAC model, activity coefficient (g residual ) is calculated from the s prole and 1/g residual is directly related to the solubility of a compound according to eqn (15). The calculated s prole of metal complexes I, II, III, IV, V and VI, shown in Fig. 5, have been compared with their experimental 46 solubility and are found to show that the metal complex having more COSMO surface is less soluble in polar aprotic solvent like acetonitrile (ACN) (i.e., II > IV > V > VI). It is also clear from our study that the metal complex with greater COSMO charge density on the metal center is found to show more solubility in the organic solvent such as acetonitrile. Our qualitative prediction of solubility in terms of the s prole of different metal complexes is consistent with the experimental solubility of those metal complexes in acetonitrile (ACN). 46 The factors affecting the solubility of a metal complex in non-aqueous media are (i) activity coefficient of a metal complex in a non-aqueous medium, (ii) molar mass, (iii) heat of fusion and, (iv) melting temperature of the metal complex according to eqn (14) and (15). According to the present study higher molar mass means lower solubility since higher molar mass is associated with the less solubility and therefore, high heat of fusion. As shown in Table 2, a metal complex with higher molar mass and bigger total COSMO surface will be less soluble in a non-aqueous medium for the same metal center with different type of electron withdrawing ligands. Also, we observed that the metal complexes with similar type of ligands and similar total COSMO surfaces have almost the same solubility in case of V(acac) 3 46 and Mn(acac) 3 ($0.60 M) 46 metal complexes, respectively. The qualitative solubility trend observed in our calculations has been explained in Fig. 6 and 7. For the metal with similar type of ligands, the results can be well explained by the ln(g 1 /g 2 ) vs. x mc , G ex (kcal mol À1 ) vs. x mc and DG mix (kcal mol À1 ) vs. x mc plot. The g 1 and g 2 are the activity coefficient of the solute and solvent molecule, respectively. The slight variation in their solubility can be easily observed in their respective plots shown in Fig. 6. The metal complex with more negative G ex (kcal mol À1 ) and DG mix (kcal mol À1 ) has more solubility in acetonitrile (ACN). To explain the solubility trend for the metal complex in presence of the electron withdrawing ligand (such as Ph), we consider the activity coefficient due to the combinatorial term to calculate the ln(g 1 / g 2 ), G ex (kcal mol À1 ) and DG mix (kcal mol À1 ). Using the information about the combinatorial term we calculated the ln(g 1 / g 2 ) vs. x mc , G ex (kcal mol À1 ) vs. x mc and DG mix (kcal mol À1 ) vs.
x mc plot for the metal complex with different type of ligands. In Fig. 7, the positive excess free energy for the metal complex V and VI correspond to the lower solubility of them in acetonitrile with respect to the IV. Therefore, our study on s prole for such metal complexes was not only able to provide a qualitative prediction of the solubility of a metal complex in acetonitrile solvent, it was able to explain the experimental solubility results too even in the absence of metal parameter in the COSMOSAC model. Also, it is noticed from our study that the effect of the entropic term due to the combinatorial interaction is signicant enough to address any trend observed in the solubility for metal complexes IV, V and VI in acetonitrile (ACN). For the metal complexes I, II and III, the residual term for the enthalpy and the combinatorial term for the entropy of the system will play a major role in order to address any trend observed in the solubility. Fig. 5, shows a comparison of the s proles of the ligands with the s proles of the metal complexes. Hence, it can be concluded that one can calculate the segment activity coef-cient for a particular metal complex even in the absence of the parameter for the metal center using COSMOSAC model.

s prole of ionic liquids and their conformational dependence
As for metal complexes, COSMOSAC model is missing optimized parameters 26 for ionic liquids. Ionic liquids have been   Fig. S-4(b). † In all cases, we found signicant differences between the s proles treated as intact ionic liquid and as solvent separated ions assuming a complete electron transferred from cation to anion shown in Fig. 8 Fig. S-6. † In all those cases, the s prole for the intact ionic liquid indicates that the separate ion model exaggerates the COSMO surface points for the anion shown in black solid line shown in Fig. S-6 in ESI. † The difference between the total s prole of separate cation and anion and the total s prole of the neutral ionic liquid is indicative towards the interaction between cation and anion of an ionic liquid shown in Fig. 8. If they did not interact with each other, the difference would be a straight line passing through zero. Also, we observed with the increasing COSMO surface of the ionic liquid the height of s prole increases for [NTf 2 ] À based ionic liquids. It is expected that the bigger molecule with small intermolecular interactions between cation and anion will be responsible for the increasing solubility of a redox active species with minimization of the Gibbs free energy of mixing which should be negative during the dissolution process. Also, treating the ionic liquid as (case I) a single molecular entity and as (case II) separate entities, it will affect the total calculated activity coefficient. In the rst case, the system will be treated as a binary mixture but in the second case it will be treated as a ternary mixture in presence of the metal complex and the calculated activity coefficient will be divided by 2 in the second case. A detailed explanation on this has been already given in ref. 38. We calculated the activity coefficient for 10 ionic liquids by considering (case I) and (case II) schemes. The results for 10 ionic liquids are given in Table 3. It is clear from our study that the different treatments of ionic liquid in the COSMO calculation will lead to different results in activity coefficient at innite  dilution. The differences are very prominent which also indicating towards the intermolecular interaction between the cation and anion of an ionic liquid. All 10 ionic liquids studied along with their COSMO surface points have been given in  ) and hence will affect their s prole, shown in Fig. 9. However, we observed signicant changes in the calculated activity coefficient at innite dilution for two different conformers because the s prole results from the COSMO surface, while the activity coefficient is a result of COSMO surface and COSMO volume. The results for activity coefficient at innite dilution for 10 ionic liquids for two different conformers are given in Table 4 and the two different conformers have been shown in Fig. S-5 in ESI. † We observed that the difference in conformers can be gauged in terms of the difference in activity coefficients. We also compare the bond energy of Conf 1 (obtained by performing ADF-COSMOSAC-2013 calculation) and Conf 2 (obtained by performing ADF-COSMOSAC-2013 calculation precede the geometry optimization and analytical frequency calculation in ADF). We observed two different conformers of same bond energy aer the COSMO calculation in homogeneous conductor while they were generated in two different ways in this study. However, to propose our model we use Conf 2 for all 10 ionic liquids for which prior geometry optimization followed by analytical frequency calculation has been done. Though we observed two different conformers of same bond energy, the calculated activity coefficient at innite dilution for two different conformers differ signicantly and this happens due to the noticeable change observed in their COSMO volume and surfaces. The results for different treatment of the ionic liquids and the different conformations have been shown in Fig. S-6 and S-7, † respectively.

Conformational dependence of s prole of metal complex
The quantum electronic structure based calculation of solubility strongly depends on the conformation of the solute and solvent species. To observe the effect in the s prole and hence on solubility of metal complexes, we have carried out our calculations for different conformations of V(acac) 3 metal complex. The details of computation have been given in Section 3. The ve different conformations of V(acac) 3 have been given in Fig. 10 and the s proles from all those conformations calculated using ADF-COSMOSAC-2013. 26 The conformational analysis using Boltzmann weight has been done to get the % populations for each conformer. We found that the % populations for Conf 5 , Conf 2 and Conf 1 were 34.71%, 33.03% and 32.27%, respectively. The details of population analysis including Boltzmann weight have been given in Table 5. For the population analysis, we consider the nal bond energy resulting from the nal COSMOSAC calculation for the ve different conformers. To see the effect of conformations on the solubility of the V(acac) 3 where, s E is the standard deviation, N is the sample size, x i is the ith value of a property and x is the mean and d E is the standard

Solubility of metal complex in ionic liquids
The calculated solubility of metal complexes V(acac) 3 and Cr(acac) 3 in ionic liquids, and correlation with the inverse of activity coefficient of metal complexes at innite dilution (1/g N ) in ionic liquids are shown in Fig. 11(a), (b), S-8(a) and (b), † respectively. From our theoretical calculations, we found that the solubility of these metal complexes will increase with increasing size of cation and anion of an ionic liquid. For all solubility calculations for V(acac) 3 (14). In the case of ideal solvation g mc ¼ 1, the percentage error in solubility is 2.12% for heat of fusion 7.17 kcal mol À1 and 4.31% for heat of fusion 5.68 kcal mol À1 calculated using eqn (14). We repeated similar calculations for Cr(acac) 3 in ionic liquids for different heat of fusions and melting temperatures already reported in the literature. [48][49][50][51][52][53] We found a similarity with the results already obtained for V(acac) 3 . The solubility calculation for Cr(acac) 3 has been given in the ESI in Fig. S-8. † The effect of different heat of fusion and melting temperature has been looked at thoroughly for V(acac) 3 3 and Cr(acac) 3 in 10 ionic liquids have been given in Table 7. The average property ( x), standard deviation (s E ) and standard error of mean (d E ) have been calculated following the eqn (16) and (17) in Section 4.
The ionic liquids containing [NTf 2 ] À anion were found to be more suitable solvents for redox active species. In the predicted solubility scale in this study, the ionic liquids containing [NTf 2 ] À show higher solubility than those having [BF 4 ] À , [OTf] À and [DCA] À anions. This can be correlated with the bigger size of the ionic liquid as seen in the corresponding s prole of those ionic liquids in Fig. 8. The bigger ionic liquid will have the maximum peak height in their s prole due to the higher COSMO surface. On the cation side, the imidazolium and pyrrolidinium cation are found as suitable ionic liquids to increase the solubility of metal complex in redox ow cell. Our theoretical predictions are in good agreement with the experimental work already reported in the work of Katayama et al. 13,[15][16][17][18] and Ejigu et al. 54 According to their studies, it is already known that the ionic liquid containing imidazolium ring and [NTf 2 ] À anion is the suitable non-aqueous solvent medium for redox active species containing transition metal center. We calculated G ex and DG mix for three different ionic liquids with solubility order:

(i) least ([Eohmim][BF 4 ]), (ii) medium ([C 8 mim][BF 4 ]) and (iii) maximum ([C 4 mim][NTf 2 ]
) to justify our solubility results of V(acac) 3 in those ionic liquids. In Fig. 12(a), the model was able to capture the asymmetric interaction present in the binary mixture solutions, the monotonic linear relation corresponds to the symmetric nature of the solution whereas the non linear interaction corresponds to the asymmetric nature of the metal complex and ionic liquid  ]) containing strong intermolecular hydrogen bond forming moieties between cation and anion shows the least metal complex solubility. This can be attributed to the high cavitation energy needed to accommodate the metal complex in this solvent and thus will reduce the solubility of metal complex in that ionic liquid and will result positive G ex and DG mix . The calculated free energy of mixing (DG mix ) and excess free energy (G ex ) for three different ionic liquids showing least, medium and maximum solubility of redox active species V(acac) 3 are shown in Fig. 12. The symmetric interactions between the V(acac) 3 4 ] are shown in the same Fig. 12(a). This is indicative to the less molecular interaction present during the dissolution process of V(acac) 3 26,55 classies the segment of the COSMO surface into three categories: (a) non hydrogen bonding, (b) hydrogen bonding from OH group and, (c) hydrogen bonding from other than OH group. A Gaussian like function has been considered to express the probability of hydrogen bonding segments where s is screening charge density and s 0 is equal to 0.007 e A À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 s prole due to the non-hydrogen 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. 26 The molecule with the ability to form intermolecular hydrogen bond between hydroxyl (OH) group of the cation and FB group of [BF 4 ] À is the least V(acac) 3 absorber because the physical absorption of V(acac) 3 molecule in those ionic liquids will result high Gibbs free energy of cavity formation in expense of breaking of intermolecular (H/F) hydrogen bonds. Such kind of (OH/OT) bond formation has been shown in the Fig. 13; plot (a) by blue line, the smaller amount of (OH/OH) has been shown by green line and the total of (OH/OT) and (OH/OH) has been shown by the red line and the total s prole p s (Total) has been shown by the black line. Also, it has been found from the study that the amount of (HB-OT) bond formation decreases from [Eohmim] [BF 4 ] to [C 4 mim][NTf 2 ] with simultaneous increase of the V(acac) 3 solubility in them. This can be correlated with the excess Gibbs free energy G ex and Gibbs free energy of mixing DG mix in Fig. 12. The molecule with the ability to form stronger hydrogen bond between the cation and anion is found to show the less solubility of V(acac) 3 in them. In the other two ILs showing medium and maximum V(acac) 3 solubility, they have no hydroxyl group (OH). Therefore all those hydrogen bond interactions in those ILs occurs due to the other than hydroxyl group (OT) interactions are weak and these results are reected in their corresponding s proles in Fig. 13 In order to investigate the effect of residual and combinatorial effect to screen the ionic liquids for a particular metal center such as V(acac) 3 , we performed our calculations for the individual effect of residual and combinatorial term on the solubility of metal complex. We observed that the effect of combinatorial term is very negligible and of the same order for all cases (1/g comb $ 1). The residual effect is sufficient enough to screen the ionic liquids for V(acac) 3 . The computed residual and Table 7 The solubility of V(acac) 3 and Cr(acac) 3 in 10 ionic liquids at different melting temperature and heat of fusion. combinatorial terms for V(acac) 3 in ionic liquids are shown in Fig. 14. Thus, based on our calculations, one will be able to do the solubility prediction of a metal complex such as V(acac) 3 just by looking at the residual interaction present in the system consisting of transition metal complex and ionic liquid using eqn (15). The 1/g Res is composed of the contribution coming from the residual and dispersion interaction in the COSMOSAC-LANL model in eqn (15) which is shown in Fig. 14. Therefore, the information about the residual interaction will be sufficient enough to screen the ionic liquids for the solubility of metal complex V(acac) 3 in them.

Dual solute effect between metal complex and ionic liquid
It has been reported that the metal complex will not behave as a sparingly soluble salt in presence of ionic liquid as supporting electrolyte and acetonitrile (ACN) as organic solvent. 1,9 The experimental studies have been done for ternary systems consist of V(acac) 3 , [Tea] [BF 4 ] and ACN and the same system with Cr(acac) 3 and Mn(acac) 3 metal complex. For all cases, it was observed that the solubility of metal complex decreases with the increase of ionic liquid concentration and this phenomena in the redox ow cell is called dual-solute effect.
We have established this effect by considering a binary system consisting of redox active species and ionic liquid for all V(acac) 3 , Cr(acac) 3 and Mn(acac) 3 metal complexes. The correlation has been established between the calculated solubility and the inverse of activity coefficient at innite dilution (1/g N ) of metal complexes with the varying mole fraction of the [Tea][BF 4 ] ionic liquid. We have found that the calculated phase diagram in absence of the third component is in good agreement with the experimental phase diagram reported at the room temperature 9 and thus was able to explain the reason behind the experimentally observed dual-solute effect. It is to be noted that the calculations have been done at room temperature to mimic the actual experimental condition used in the article. 1 From Fig. 15(a), (b), (d), (e), (g) and (h) of all species in the binary mixture, it is found that the  dual-solute effect is a consequence of Gibbs-Duhem relationship 28 between the two species in a binary mixture in the COSMOSAC model. The solubility phase diagrams already found in the work of ref. 9, has been shown in Fig. 15(a), (b), (d), (e), (g) and (h) are in good agreement with the experimental 9 phase diagram reported for V(acac) 3 , Cr(acac) 3 and Mn(acac) 3  For the least and medium absorber ionic liquids, we observed that the activity coefficient of the solute is increasing with the increase of the mole fraction of the ionic liquid while the opposite behavior has been observed for the maximum absorber ionic liquid which can be explained in terms of the symmetric and asymmetric interaction ( Fig. 16(a)) between the solute and solvent species and can be directly correlated with the Gibbs-Duhem relationship and solubility ( Fig. 16(b), (c), (d) and (e)) observed in each case. Also, these theoretical ndings are conrmed by the calculations of ln g 1 /g 2 , G ex and, DG mix as a function of mole fraction of V(acac) 3 in Fig. 12. Fig. 15 The calculated phase diagrams have been shown in (a) and (b) for V(acac) 3 , in (d) and (e) for Cr(acac) 3 and in (g) and (h) for Mn(acac) 3 metal complexes. In (c), (f) and (i), the Gibbs-Duhem relationship for the three binary systems have been shown.

Conclusion
We have used COSMOSAC-LANL model to predict solubility of metal complexes in ionic liquids and to screen the ionic liquids for a particular metal complex where the metal complex has been treated as a solute and the ionic liquid as a solvent in a binary mixture. Also, we have shown that this model can do qualitative prediction of solubility of different metal complexes with same metal center but with different type of ligands for a particular solvent by comparing their s proles using the model. The proposed new model was able to explain the similar experimental solubility results observed in the case of V(acac) 3 , Cr(acac) 3 and Mn(acac) 3 in acetonitrile solvent. The bigger size of ionic liquids and metal complexes have been correlated with the solubility of metal complex to get a qualitative idea on solubility trend. We have applied our theoretical model on a series of chromium complexes to get qualitative prediction of their solubility in acetonitrile solvent using their s prole information because this information is used for the calculation of solubility of the metal complex in the ionic liquid. We found our theoretical results are in good agreement with the experimental results. 46 According to our calculations, the ionic liquids containing imidazolium cation and [NTf 2 ] À are found to show greater solubility of redox active species in them and our results are found to follow the trend already reported in the experimental paper. 13,[15][16][17][18]54 We calculated the excess free energy (G ex ) and free energy of mixing (DG mix ) for least, medium and maximum absorber of redox active species in ionic liquids and found that the free energy of mixing is less negative for [ ]. The present model was also able to explain the solvation mechanism in those systems. This model can be used to screen the suitable ionic liquid for a particular metal complex and metal complexes with different ligands and with same metal center against a particular solvent to increase the ionic conductivity, solubility and thus energy density of a redox ow cell where a compatible pairing between metal complex, organic solvent and ionic liquid is an important factor.

Conflicts of interest
There are no conicts to declare.

Appendix A
At innite dilute condition, when x 1 / 0 and x 2 / 1, the eqn (6) will be reduced to ln g N 1 ðasymÞ/ a 1 þ b 1 RT for solute species (1). Now aer adding and substituting the value of (a 1 + b 1 ) in the expression of ln g N 1 , one will get ln g N 1 ðasymÞ ¼ A 12 RT and hence, Similarly, when x 2 / 0 and x 1 / 1, the eqn (7) will be reduced to ln g N 2 ðasymÞ ¼ A 21 RT for the solvent species (2) and therefore, Therefore, at innite dilution, one can calculate Margules parameters from the activity coefficients at innite dilution from the above two equations. Since, in this article the Margules parameters are called from the activity coefficient at innite dilution calculated using COSMOSAC-2013 model, therefore one can write and ln g N 2 ðCOSMOSACÞ ¼ Substituting the value of A 12 and A 21 in eqn (19) and (20), one will get, and ln g N 2 (asym) ¼ ln g N 2 (COSMOSAC).
Therefore, subtracting eqn (25) from eqn (28), one will get the relationship between the two models which is valid at innite dilution and that is ln g(COSMOSAC-LANL) À ln g(COSMOSAC-2013) ¼ ln g(comb). (29)