Computation of gas solubilities in choline chloride urea and choline chloride ethylene glycol deep eutectic solvents using Monte Carlo simulations

Deepeutecticsolvents(DESs)areconsideredasgreenalternativestoroomtemperatureionicliquids(RTILs),due to their lower-cost synthesis and more environmentally friendly nature. In this work, Monte Carlo (MC) simula-tionshavebeenusedtocomputethesolubilitiesofCO 2 ,H 2 S,CH 4 ,CO,H 2 ,andN 2 incholinechlorideurea(ChClU) and choline chloride ethylene glycol (ChClEg) DESs. Due to the strong intermolecular interactions of DESs, lead-ingtohighviscosities,MCsimulationspresentsigni ﬁ cantchallengeswithrespecttosystemequilibrationandsol-ute molecule insertions. The Continuous Fractional Component Monte Carlo (CFCMC) method has been used with our open-source code, Brick-CFCMC, to improve molecule insertions and equilibration of the system, and directly compute the excess chemical potential and solubility (in terms of the Henry constant) of the gas mole- cules in the DESs. Pure DES properties, such as density and radial distribution functions (RDFs), were well reproduced by MC simulations. The solubilities of gases were, however, underestimated by the CFCMC simula- tions compared to available experimental data from literature. The order of solubilities of the different gases in ChClU at 328 K was obtained as H 2 S > CO 2 > CH 4 > H 2 > CO > N 2 , which reasonably agrees with experimental data from literature. The OPLS force ﬁ eld resulted in larger average Henry constants in ChClEg, compared to the GAFF force ﬁ eld, implying the better suitability of the GAFF force ﬁ eld for the calculations. Smaller ionic charge scaling factors were shown to increase the solubilities of the gases in the DESs, but result in lower densities. The differences between the computed Henry constants from MC simulations and experimental data from liter- aturemaybecausedbytheunsuitabilityoftheusedforce ﬁ eldparametersoftheDESsincombinationwiththose of the solute gases. Nonetheless, experimental data from literature is scarce (except for CO 2 ) and in some cases contradictory, which makes the comparison with the computational results dif ﬁ cult.


Introduction
In recent years, deep eutectic solvents (DESs) have received a growing attention as 'green' alternatives to the conventional volatile organic compounds (VOCs), due to superior properties, such as low vapor pressure and non-flammability [1][2][3][4][5][6]. DESs and room temperature ionic liquids (RTILs) share many physical and thermal properties, such as high solvation with respect to many solutes, low vapor pressure, wide temperature range for liquidity, and low melting point [5,[7][8][9][10]. These properties make DESs suitable media for separation processes. Easy and inexpensive preparation with no further purification requirement, non-toxicity and biodegradability are some of the potential benefits that make DESs more attractive solvents compared to many RTILs [1,11,12]. DESs are most commonly synthesized by mixing a hydrogen bond donor (HBD) with a hydrogen bond acceptor (HBA), forming a eutectic mixture with a melting temperature significantly lower than those of the individual components [8]. This depression in melting point has been attributed to the extensive hydrogen bond network formed between various moieties in the mixture [3,9]. The molar ratio in which the HBA and HBD components are mixed is referred to as the eutectic molar ratio and plays a crucial role in determining the physico-chemical properties of DESs [1]. In many well-known hydrophilic DESs, the HBA component is an ion pair, such as choline chloride. However, DESs can, in principle, form by mixing non-ionic species [10,13,14]. More details on DESs and their potential applications can be found in literature [1,[8][9][10][15][16][17][18].
In this work, MC simulations were used to obtain the solubility of CO 2 , H 2 S, CH 4 , CO, H 2 , and N 2 molecules in choline chloride urea (ChClU) and choline chloride ethylene glycol (ChClEg) DESs with HBA:HBD molar ratios of 1:2 (corresponding to the eutectic point), as well as the densities and RDFs of the pure DESs. In these DESs, choline chloride is the HBA component, while urea and ethylene glycol are the HBD components of ChClU and ChClEg, respectively. The gas solubilities were computed in terms of the Henry constant at infinite dilution. As several of the aforementioned gases are not very soluble in DESs, solubility measurements for these gases can be challenging. The solubility measurement data for DESs in literature are mostly limited to CO 2 and SO 2 . The reported experimental solubilities from different sources are in some cases inconsistent. For instance, the mole fraction based Henry constant of CO 2 in ChClU has been reported as ca. 57 MPa (5 MPa molality based) at 328 K by Mirza et al. [63], while several other studies have published values in a range of 16 MPa to 20 MPa [25,60,65,74]. This inconsistency may be due to differences in experimental methods and conditions used by the various research groups. Notably, a variety of pressure ranges are used for these solubility measurements, e.g., 8.5 bar to 125 bar [74], 6 bar to 45 bar [65], 0.1 bar to 2 bar [60], 3 bar to 60 bar [25], and 0.4 bar to 1.5 bar [63]. Therefore, it is possible that in some of the high pressure measurements of the absorption isotherm, the Henry regime was not reached. Furthermore, the water content of the DES varies between experimental studies. The presence of water in DESs may affect the intermolecular interactions and physical properties of DESs [75][76][77]. Studies have suggested that the water content of a DES may behave as an antisolvent and adversely influence the CO 2 solubility [78,79]. Another example of the inconsistent solubilities in literature is the Henry constant of CH 4 in ChClU at 328 K, for which Liu et al. [60] reported a mole fraction based value of 231 MPa (extrapolated value), while Xie et al. [65] obtained a value of 47 MPa. It is also possible that rather than differences in measured solubilities, the inconsistency lies in the reporting of the data, e.g., by assuming different units and definitions, without clearly stating these assumptions. Due to the limited availability of experimental solubility data in literature for DESs, and the inconsistency of the data, molecular simulation may be considered as an alternative for these experiments.
This manuscript is organized as follows: The next section is devoted to the description of the simulations details. This consists of the force field parameters used for the various molecules and the computational methods to obtain densities and Henry constants. Subsequently, the results from the simulations are presented, discussed and compared with the limited experimental data from literature. Finally, conclusions are provided regarding MC simulations of DESs for solubility calculations.

Force field parameters
The non-polarizable force field parameters by Perkins et al. [37,38], based on the Generalized AMBER Force Field (GAFF) [80] were used for ChClU and ChClEg. The force field parameters consist of bonded terms (bond stretching, angle bending and torsion), and Lennard-Jones (LJ) and electrostatic non-bonded terms. Following the approach of Perkins et al. [37,38], ionic charges were scaled by 0.8 and 0.9 for ChClU and ChClEg, respectively, to take effective polarization into account. The ionic charge scaling approach is often used in IL simulations, as it yields a better agreement between the simulation results and experimental data [81]. The 1-4 intramolecular interaction energies were scaled according to the AMBER [82] force field, i.e., scaling factors of 0.5 for LJ and 0.833 for electrostatics. To examine the effect of force field on the computed solubilities, the non-polarizable force field parameters by Doherty and Acevedo [39], based on the OPLS force field [83], were also used for ChClEg. This force field also comprises the aforementioned non-bonded and bonded terms. For the OPLS model, the ionic charges were scaled by 0.8, according to the original publication [39], and the 1-4 LJ and electrostatic intramolecular energies were scaled by 0.5 [83]. For both force fields, LJ parameters were set to ε/ k B =0.5K and σ=0.1 Å for the unprotected hydrogen atoms in hydroxyl groups to prevent atomic overlaps [84]. All bonds were kept fixed at the equilibrium bond lengths during simulations. The improper torsion potentials were not taken into account. MD test simulations showed that this had little effect on the density of ChClU using the GAFF force field. From these test simulations, the density at 328 K was computed as 1.199 kg m −3 and 1.217 kg m −3 , with and without the improper dihedrals, respectively, showing a relative difference of ca. 1.5%. All solute molecules were modeled as rigid objects. The Transferable Potentials for Phase Equilibria (TraPPE) force field was used for CO 2 , CH 4 (united atom), N 2 , and H 2 S [85,86]. The two-site force field parameters by Cracknell [87] were used for H 2 , and the three site model by Martín-Calvo et al. [88] was used for CO. Force field parameters for all the molecules are listed in the Supporting Information (SI). Long-range electrostatic interactions were computed with the Ewald summation method [89,90]. The Ewald parameters were set based on a relative precision of 1 × 10 −6 . A cutoff radius of 10.0 Å was used for both LJ and shortrange electrostatic interactions. For ChClU, the LJ potential was shifted at the cutoff radius, and no analytic tail corrections [91] were used, whereas for ChClEg, analytic tail corrections were used, but the LJ potential was not shifted. The choices regarding the shifting of the LJ potential and the use of analytic tail corrections were made to obtain a closer agreement of the computed densities of the pure DESs with experimental densities reported in literature [21,[92][93][94]. To accurately compute the Henry constants, however, analytic tail corrections were used for solute molecules in both DESs to account for the long-range solvent-solute LJ interactions. Lorentz-Berthelot and Jorgensen mixing rules were applied to calculate the LJ interactions between non-identical atomtypes for the GAFF and OPLS models, respectively.

Computational methods
The densities of ChClU and ChClEg, as well as solubilities, i.e., Henry constants, of the gas molecules in these DESs were calculated from MC simulations in the isobaric-isothermal (NPT) ensemble. The densities were computed at 1 bar and various temperatures (308 K-338 K). RDFs were computed between various atoms in the system during the same simulations as for density computation. Solubilities were calculated in a separate set of simulations at a temperature of 328 K and a hydrostatic pressure of 1 bar. 50 HBA and 100 HBD molecules were used for both DESs, corresponding to the HBA:HBD eutectic molar ratio of 1:2. Initial configurations were generated at a lower density and the simulation box was compressed to the equilibrium density at the specified pressure. The equilibrium box size in these simulations was typically around 27 Å. After the system was equilibrated, the production cycles followed for computing average properties. During each MC cycle, trial moves were performed with fixed probabilities. The number of trial moves per MC cycle was equal to the total number of molecules in the system. For each density data point, 10 independent simulations were run for 1 × 10 6 equilibration and 1 × 10 6 production cycles. Solubility simulations comprised 20 independent runs, for which up to 4 × 10 6 equilibration and 12 × 10 6 production cycles were used. For the calculation of solubilities, block averaging was performed, where the results of the production runs were divided into 4-7 blocks, depending on the system. For each block, the average Henry constant was computed over the 20 independent runs. The overall mean and standard deviation values were then computed for these block averages. Trial moves for thermalization of the system included translation and rotation of the molecules, volume change of the simulation box and changing the internal configuration of the molecules. The internal configuration of molecules was altered using random changes in bond angles and dihedrals of the molecules. Maximum displacements for all thermalization trial moves were adjusted to allow for 50% acceptance probabilities.
All simulations were carried out using the open-source MC code, Brick-CFCMC [95], which is developed in our group for phase and reaction equilibria, and is particularly optimized for high-density liquidphase simulations with difficult equilibration, as in the case of ILs and DESs. Brick-CFCMC applies a modified version of the CFCMC method [96,97] by Poursaeidesfahani et al. [98] that enables direct computation of excess chemical potentials (and Henry constants) for various molecules in the system. The CFCMC method was not used during the simulations for density (and RDF) computation and was only applied in the solubility simulations. In this method, one or more extra molecules, referred to as fractional molecules, are introduced to the system. All other molecules are referred to as whole molecules. The interactions of a fractional molecule with other molecules are scaled by a coupling parameter, λ, which can have a value in the range [0,1]; a value of 0 means that the fractional molecule has no interactions with other molecules (ideal gas molecule) and a value of 1 indicates that the fractional molecule is fully interacting (whole molecule). The value of λ is changed in the course of the simulation by Δλ, essentially inflating or deflating the fractional molecule. The value of Δλ is taken randomly from a uniform distribution in the range [−Δλ max , Δλ max ], where Δλ max is fixed (here at 0.2) during the simulation. Two additional trial moves were performed in the simulations, next to thermalization and λ-change trial moves: 1. The fractional molecule was reinserted at a random position in the box, without changing the value of λ. 2. An identity change trial move was performed, where the fractional molecule was turned into a whole molecule and a randomly selected whole molecule of the same type was changed to a fractional molecule with the same previous value of λ.
The observed probability distribution of λ is generally non-uniform due to the free energy barriers encountered when the value of λ changes. Therefore, a biasing weight function, W(λ), is constructed to overcome these free energy barriers and sample the λ-space with equal probability. The Boltzmann average of any observable A is then computed using [98]: The Wang-Landau algorithm [99,100] was used in the equilibration runs to construct an initial biasing weight function. The iterative scheme [101] was applied in between consecutive production runs to further modify the weight function and yield a more uniform observed probability distribution of λ. 1000 bins were used to obtain a histogram of λ values and thus the probability of occurrence for each value. Two fractional molecules were used in each solubility simulation, one of the HBD and one of the solute gas. The fractional molecule of the HBD was used to enhance the equilibration of the system (by λ-change and identity change trial moves), while the fractional molecule of the solute was used to compute the solubility of the solute in the DES. To increase the interactions of a fractional molecule from λ=0 to λ=1, first, the LJ interactions were switched on until λ = λ switch (where λ switch was set to 0.8), by using the scaling parameter λ LJ . In this range of λ, electrostatic interactions of the fractional molecule were switched off. In the range λ = λ switch to λ=1, the electrostatic interactions were switched on by using the scaling parameter λ el , while the LJ interactions remained fully switched on. A similar procedure was used to reduce the interactions of the fractional molecule from λ=1 to λ=0, where the electrostatic interactions were scaled in the range [λ switch , 1] (while λ LJ = 0), and the LJ interactions were scaled in the range [0, λ switch ] (while λ el = 1). The intermolecular LJ energy, U LJ , between interaction sites i and j was scaled according to [102,103]: where r ij is the distance between i and j, σ ij and ε ij are the LJ parameters, and α, b and c are the LJ scaling parameters. These scaling parameters are often set to α=0.5, b=2, and c=6 in CFCMC simulations [98,[104][105][106]. In this work, however, the values for these parameters were set to α=0.0025, b=1, and c=48. This set of parameter values have been reported to minimize the statistical variance of the derivative of the total energy with respect to λ, resulting in a more efficient sampling of λspace [103]. Similar to the approach of Hens and Vlugt [105], the realspace and exclusion terms of the electrostatic energy were scaled by multiplying each of the charges in the fractional molecules by λ el , and replacing r ij by r ij + A(1 − λ el ) (where A was set to 0.01 Å) to avoid singularities at small values of r ij . Only the intermolecular LJ and electrostatic interactions were scaled in the simulations and no intramolecular interactions (bond bending, torsion, LJ, and electrostatic energies) were scaled.
The molality based Henry constant of solute i, H m , is defined as [72,107,108]: in which f i is the partial fugacity of the solute molecule in the gas phase, m i is the molality of the solute in the solution, in mol kg −1 , and m 0 is set to 1 mol per kg of the solvent. The unit of the Henry constant in Eq. (3) is therefore the same as the unit for fugacity, e.g., MPa. The fugacity of solute i is the product of its fugacity coefficient, ϕ i , and the total pressure, P.
At low pressures, such as the pressures used in this work, ϕ i → 1. Therefore, f i can be replaced by the partial pressure, P i , in Eq. (3). Moreover, due to the negligible vapor pressure of DESs, the vapor phase may be assumed to contain only the solute gas molecules. P i is thus equal to the total pressure (P), imposed in the simulations. The Henry constant defined in Eq. (3) can be computed as a function of the excess chemical potential of the solute at infinite dilution, μ ex, ∞ , as [89,109]: where R is the universal gas constant (in Jmol −1 K −1 ), T the temperature (in K), and ρ the mass density of the pure DES (in kg m −3 ). The resulting Henry constant is thus in units of Pa. It is also possible to define the Henry constant in terms of mole fraction, where m i /m 0 is replaced by the mole fraction of the solute, x i , in Eq. (3), and the right-hand side of Eq. (4) is divided by the molar mass of the DES (in kgmol −1 ). To compute the Henry constant from Eq. (4), the excess chemical potential of each solute was obtained at infinite dilution from the unbiased (Boltzmann sampled) probability distribution of λ of the fractional molecule of the solute, when λ approaches 0 (p(λ → 0)) and 1 (p(λ → 1)) [98]: To verify the Henry constant calculations, the molality based absorption isotherm of CO 2 in ChClEg (using the GAFF model) was computed at 328 K from CFCMC simulations in the expanded osmotic ensemble [52,53,96,110]. In the osmotic ensemble, the system is considered in equilibrium with an ideal gas reservoir. The temperature and hydrostatic pressure (of the system and the reservoir), as well as the chemical potential of the solute in the reservoir and one extensive variable of the system, here the number of DES molecules, are kept fixed. The number of solute gas molecules in the mixture and the volume of the system are changed to ensure phase equilibrium between the system and the reservoir. A pressure range of 0.1 bar to 10 bar was used to compute the CO 2 loading. The reservoir was assumed to be composed only of the solute gas (with a mole fraction of 1), as DESs have negligible vapor pressure. The chemical potential of the solute at each pressure was specified in terms of its fugacity coefficient, which was determined from the NIST REFPROP [111] software, based on the Peng-Robinson EOS [112]. A fractional molecule (osmotic fractional molecule) of the solute gas was used for molecule insertions in (and deletions from) the mixture. An additional fractional molecule (NPT fractional molecule) of the solute gas was used in the system to allow for the identity change trial move and thus enhanced system equilibration. The identity change and reinsertion trial moves performed in the osmotic ensemble were similar to the NPT ensemble trial moves described earlier in this section. The λ-change trial moves were performed in the osmotic ensemble in the same way as in the NPT ensemble for when 0<λ<1. However, if after a λ-change trial move the value of λ for the osmotic fractional molecule became larger than 1, the fractional molecule was transformed into a whole molecule and a new fractional molecule was randomly inserted in the simulation box with a scaling parameter of λ − 1 [96]. Similarly, if the value of λ dropped below 0, the fractional molecule was deleted from the simulation box and a randomly selected whole molecule was turned into a fractional molecule with a value of λ + 1 for the scaling parameter [96]. In the latter case, if no whole molecules were present in the simulation box, the trial move was rejected. For more details about the MC simulations, the reader is referred to Ref.

Density and radial distribution functions
The densities of ChClEg and ChClU were computed at various temperatures without the use of the CFCMC method. The resulting densities are shown in Fig. 1 as a function of temperature for the two DESs. The results are compared with the densities obtained from MD simulations of Perkins et al. [37,38] with the same force fields, as well as with experimental data [21,[92][93][94]. It can be observed that the results from MC slightly deviate from the MD simulations of Perkins et al. In the case of ChClEg, the densities from MC simulations are slightly larger than the values by Perkins et al., and for ChClU, the densities from MC simulations are slightly smaller. As bond stretching and improper torsion are not yet implemented in Brick-CFCMC, these potentials were neglected in the MC simulations. The effects of the bond stretching and improper torsion exclusions, combined with the use of tail corrections, shifting of the LJ potential, and dissimilar cutoff radius and system size may be the cause of the slight differences in density between the MC simulations and the results of Perkins [37,38] (red circles), and experimental data of Yadav and Pandey [92] (blue squares), Leron and Li [93] (green diamonds), Yadav et al. [21] (blue squares), and Leron et al. [94] (green diamonds). The small differences observed between the densities obtained from MC in this work and the MD simulations of Perkins et al. [37,38] may be due to the exclusion of improper torsion and bond stretching potentials from the MC simulations, the use of tail corrections, shifting of the LJ potential, and difference in the cut off radius and system size. similar experimental densities for the DESs [79,113,114]. This indicates that the densities of DESs can be accurately computed using MC simulations, even when more advanced techniques such as CFCMC are not used.
RDFs between various atoms were computed for both DESs with the GAFF force field parameters (Fig. 2). The results agree with the RDFs reported in literature using MD simulations with the same force fields [37,38,115]. For instance, for the HO choline − Cl RDF in ChClEg, both the MC simulations in this work and the MD simulations of Perkins et al. [38] show a first peak at ca. 2.3 Å with an intensity of ca. 11. Therefore, MC simulations are able to reproduce, as accurately as MD simulations, the liquid structure of the DESs, without the need to use the CFCMC technique. Considering the relatively high viscosity of the DESs at the simulation temperature of 328 K (24 cP for ChClEg, and 95 cP for ChClU [22]), this is an important finding. However, it is expected for such calculations to be more challenging at lower temperatures, closer to the glass transition temperature, where the viscosities of DESs are significantly higher [15,92,116]. The RDFs computed from the MC simulations of ChClU are in agreement with the model fitted experimental data (Neutron diffraction) of Hammond et al. [117]. For instance, two peaks are obtained from the MC simulations for the O urea − Cl RDF, at 4.5 Å with an intensity of 1.5 Å, and at 5.2 Å with an intensity of 2.2 Å. Hammond et al. [117] reported comparable peak intensities, although the position of the first peak is reported slightly lower (ca. 3.8 Å) compared to the MC simulations. Several other computational studies are available that report RDFs for ChClU and ChClEg, based on MD and density functional theory (DFT) simulations [33,47,49,50]. Differences can be observed to various extents between the MC computed RDFs and the RDFs reported in literature, due to the dissimilarity in modeling methods and/or parameters (e.g., force fields). For instance, in the case of ChClU, the MC computed RDFs are in agreement with the MD results of Kumari et al. [50], using the CHARMM 36 force field [118,119], except for some differences in the first peak intensities. The first peak intensities for the N urea − Cl, H urea − Cl, and HO choline − Cl RDFs, for instance, were computed from MC simulations as ca. 4.9, 4, and 5.2, whereas Kumari et al. reported the first peak intensities of these RDFs as ca. 6.8, 6, and 11. The position of peaks and the intensities of the other peaks (higher solvation shells) are, however, in agreement between the MC simulations and the results of Kumari et al. for these RDFs. Differences between RDFs are more prominent when comparing the MC simulation results with the MD results reported by Sun et al. [49], where force field parameters based on OPLS/ AMBER [120] and an ionic charge scaling factor of 0.9 were used. For ChClEg, the RDFs from MC simulations in most cases differ in the intensity of the first peak from the RDFs of Ferreira et al. [33] using MD simulations with force field parameters based on OPLS [83,121,122], and an ionic charge scaling factor of 0.8. For instance, the first peak intensity of the HO eg − Cl RDF was computed as 11.4 from MC simulations, while Ferreira et al. reported a value of ca. 7. It is noteworthy that the ionic charge scaling factors in the aforementioned studies were adjusted for the best agreement of the simulation results with experimental data. Therefore, depending on the accuracy of the force field used in the simulations, different scaling factors were used by different authors. The differences between the RDFs from various studies highlight the fact that accurate modeling of the liquid structure of DESs using MC/MD simulations significantly depends on the use of optimal force field parameters.

Henry constants
The Henry constants of CO 2 , H 2 S, CH 4 , CO, H 2 , and N 2 in ChClU and ChClEg DESs were computed using the CFCMC method in the NPT ensemble. Equilibration of the simulations was verified by monitoring the changes in density and total energy of the mixture. The liquid structure of the mixture was investigated by computing RDFs. The obtained RDFs were in agreement with the RDFs computed in the density simulations (no CFCMC). This means that the additional solute and HBD fractional molecules in the system did not change the structure of the DES. This is consistent with the observations of Rahbari et al. [106] that the presence of a small number of fractional molecules (up to ca. 1% of the total number of molecules) does not change the thermodynamic properties of the system. The average observed probability distribution of λ was investigated for the solute fractional molecules at the end of each production run as a measure of λ-space sampling. Flat observed probability distributions of λ were obtained for all the solutes in the DESs (not shown here), indicating that the sampling of λ-space was sufficiently performed. As explained in Section 2.2, to obtain flat λ probability distributions (and reduce uncertainties in the Henry constants), long equilibration and production runs (several million MC cycles) were required, owing to the high viscosity and density of the DESs. The obtained weight functions for all the solute fractional molecules in the ChClU mixture, averaged over all independent runs, are shown in Fig.  3. The difference in weight function between λ=0 and λ=1 represents the free energy change for insertion/deletion of the solute. By examining the weight functions of the different solutes in ChClU, in Fig. 3, it becomes apparent that CO 2 and H 2 S are much more soluble than CH 4 , CO, H 2 , and N 2 . It can be observed that for the insoluble gases CH 4 , CO, H 2 , and N 2 , the free energy barrier of insertion is larger in the LJ scaling range (0 to 0.8) compared to the electrostatic scaling range (0.8 to 1). This free energy barrier is encountered when the interactions of the fractional solute molecule with its surrounding molecules are gradually scaled by λ, i.e., the molecule is gradually inserted (or deleted). The free energy barrier of insertion/deletion thus depends on the λ-scaling pathway used in the simulations, for instance by setting different values for the constants α, b, and c, in Eq. (2). The models used for CH 4 and H 2 have no partial charges and therefore, have flat weight functions in the electrostatic scaling range. The LJ interactions are also observed to impose a larger free energy barrier for the insertion/deletion of CO 2 and H 2 S compared to the electrostatic interactions, although the differences are smaller than in the case of insoluble gases. Similar λ probability distributions and weight functions were obtained for the dissolution of gases in ChClEg (not shown here).
The computed molality based Henry constants from the CFCMC simulations are listed in Table 1 for both DESs. For ChClEg, the obtained values from both GAFF and OPLS force fields are provided in the table. Except for CO 2 , experimental data are very scarce for the solubilities of the studied gases in ChClU. To the best of our knowledge, no experimental data exist for the solubility of these gases (except for CO 2 ) in ChClEg. Some of the experimental data found in literature are also listed in Table  1. In case the reported Henry constants were on a mole fraction basis, a conversion to molality based values was performed, using the molar mass of the DESs (86.57 g mol −1 for ChClU, and 87.91 g mol −1 for ChClEg). It can be observed that for the Henry constant of CO 2 in ChClU, the computed value from CFCMC (6 MPa) is a factor of 4 larger than the experimental data by Xie et al. [65] and Liu et al. [60]. Leron et al. [25] and Li et al. [74] have also reported the Henry constant of CO 2 in ChClU to be in the range 1.4 MPa to 1.7 MPa (molality based). However, as shown in Table 1, the Henry constant by Mirza et al. [63] is much larger than the other experimental data, and is closer to the computed value from CFCMC. As discussed in Section 1, such inconsistencies may be due to the different pressure ranges used for the solubility measurements. Furthermore, the water content of the DES used by Mirza et al. [63] (wt% 2.4) is higher than in the other studies, which could be the reason for the larger reported Henry constant. However, at a similar water content of ChClU (wt% 1.9), Xie et al. [79] reported a Henry constant of 1.75 MPa (molality based) for CO 2 . The authors also showed that upon an increase in the moisture content of ChClU from wt% 1.9 to wt% 9.1, the Henry constant increased to 1.93 MPa, which is still lower that the value by Mirza et al. [63]. This contradiction, in addition to the scarcity of data for other gases, makes the comparison between simulation results in this work and experimental data difficult. The computed Henry constants for other gases are also significantly larger than the experimental data (up to a factor of 88 for N 2 ). As can be observed in Table 1, the reported experimental Henry constants for CH 4 differ by a factor of 5. This is an indication that more experiments are required to confirm the solubility values for such insoluble gases in ChClU. Nonetheless, it can be observed that the relative order of the Henry constants for the various solutes agrees reasonably well with the experimental data. The solubilities of gases from simulations are in the order H 2 S > CO 2 > CH 4 > H 2 > CO > N 2 , whereas experimental data indicate the solubilities to be H 2 S > CO 2 > CH 4 > N 2 > H 2 > CO at the same temperature, differing only in the relative solubility of N 2 . It is important to note that the relative magnitude of the Henry constant for the various gases is consistent with the observations regarding the weight functions of the gases (Fig.  3), as a representation of the free energy change of dissolution. Using the GAFF force field parameters, the computed Henry constant of CO 2 in ChClEg is equal (within the uncertainties) to the Henry constant in ChClU. The Henry constant of CO 2 in ChClEg is nonetheless in better agreement with the experimental data. It can be observed that the average Henry constants obtained using the OPLS parameters are, for all the gases except H 2 S, larger than the values computed using the GAFF parameters. The difference between the average Henry constant of CO 2 in ChClEg from simulations and the reported experimental data is thus larger when the OPLS parameters are used, implying better suitability of the GAFF force field for the calculations. For H 2 S, the GAFF and OPLS force fields result in comparable Henry constants, although the average Henry constant is slightly larger when the GAFF force field is used. The difference in the computed Henry constants using the OPLS and GAFF force fields is larger for the more insoluble gases. To draw final conclusions on which force field predicts the solubilities more accurately, experimental data are required for the solubility of these gases in ChClEg. The solubilities of gases in ChClEg, using both GAFF and OPLS force fields, have the same relative order as in ChClU, with H 2 S as the most soluble gas, and N 2 as the least soluble gas. The computed Henry constants of gases in ChClEg are to some extent smaller than in ChClU, when the GAFF force field is used for both DESs, suggesting ChClEg to be a better solvent for these gases. It is also noticeable that the uncertainties in the computed Henry constants (values in parentheses in  Table 1 Computed average Henry constants (molality based) at 328 K for various gases in ChClU and ChClEg DESs from CFCMC simulations in the NPT ensemble, using the GAFF [37,38] and OPLS [39] (only for ChClEg) force fields. The standard deviations are indicated in parentheses. Experimental data available in literature [60,63,65,72] Table 1) are smaller for ChClEg than for ChClU, using the GAFF force field for the DESs, possibly due to the weaker intermolecular interactions of ChClEg and enhanced sampling of chemical potentials of the solutes. It can also be observed that the uncertainties in the Henry constants increase for less soluble molecules, as it becomes increasingly difficult to insert such molecules in the system. When the OPLS force field is used for ChClEg, however, the computed average Henry constants are generally larger than for ChClU (except for H 2 S), implying ChClU to be a superior solvent for the gases. It would be interesting to see, with future experimental data, which of the predictions by the two force fields will be more accurate.
To further examine the methodology used for computing the Henry constants of gases in the NPT ensemble, solubility simulations were carried out in the osmotic ensemble for CO 2 in ChClEg, using the GAFF force field for the DES. Thereby, the absorption isotherm of CO 2 was computed at 328.15 K as a function of the hydrostatic pressure. The resulting isotherm is presented in terms of dimensionless molality (m/m 0 ) in Fig.  4. As can be observed in the figure, the solubility of CO 2 increases almost linearly when the pressure is elevated. The inverse of the slope of the isotherm at infinite dilution (m/m 0 → 0) provides the Henry constant. The Henry constant was obtained by a linear fit to the low pressure part of the isotherm (i.e., first five points). The fitting (with R 2 =0.999) resulted in a Henry constant of 5.5 MPa, which is in agreement with the value obtained from the NPT ensemble simulations. The lines corresponding to the Henry constants obtained from both NPT and osmotic ensemble simulations are presented in Fig. 4. The difference between the two lines is small when m/m 0 → 0, and increases when the pressure/molality increases. This is due to the fact that the Henry constant from NPT simulations was obtained at infinite dilution in accordance with its definition given by Eq. (3). Experimental solubility data by Leron and Li [72], and Mirza et al. [63] are also presented in Fig. 4. The experimental data indicate a higher CO 2 loading and a lower Henry constant (as implied by the inverse of the isotherm slope) in ChClEg, compared to the calculated values from the osmotic ensemble simulations. This is consistent with the results obtained from the NPT simulations ( Table 1). The solubility data of CO 2 in the commonly studied [bmim] [Tf 2 N] IL [123] are also shown in Fig. 4 for comparison. It can be observed that the experimental solubility of CO 2 in [bmim][Tf 2 N] is larger than the experimental value in ChClEg. This is also qualitatively consistent with the simulation results.
Overall, there is a significant difference between the computed Henry constants of gases in the DESs (for ChClU in particular) and experimental data in literature. Aside from the uncertainties revolving the experimental data, it is possible that the force field parameters used for the DESs, although yielding accurate thermodynamic properties for pure DESs [37][38][39]115], are not optimal for use in combination with the models used for the solute gases. A particularly influential parameter in the IL and DES force fields is the ionic charge scaling factor, which may significantly affect the computation of thermo-physical and structural properties [36,37,124,125]. A smaller value for the scaling factor (farther from 1) leads to weaker ionic interactions. Weaker ionic interactions result in lower densities of the IL/DES and thus larger free volume required for the dissolution of solutes. Therefore, it is expected that reducing the ionic charge scaling factor increases the solubilities of gases in the solvent (decreases the Henry constants). Conversely, a larger value for the scaling factor (closer to 1) increases the strength of cohesive ionic interactions and therefore yields higher densities and lower gas solubilities. This scaling factor depends on the specific force field parameters used in the simulations and is often tuned to obtain better agreement of the simulation results with experiments for pure DESs/ILs [37,39,124,126], in contrast to DES/IL mixtures with other compounds. It is thus possible that the scaling values used for the pure DESs are too large (and the cohesive interactions too strong) to allow for the dissolution of solutes. To investigate the influence of the ionic charge scaling factor, the Henry constants of CO 2 , H 2 S, and CH 4 in ChClU and ChClEg were computed in the NPT ensemble, using the GAFF force field for the DESs and smaller values for the ionic charge scaling factors. The charge scaling factor was reduced from 0.8 to 0.7 for ChClU, and from 0.9 to 0.8 for ChClEg. The results are presented in Fig.  5, where the Henry constants with the original charge scaling factors are also shown for comparison. It can be observed that, as expected, the reduction of the ionic charge scaling factor results in a decrease in the Henry constants of all the solutes (increased solubility). The Henry constants in ChClU with the charge scaling factor of 0.7 are obtained as 4.1 MPa, 1.7 MPa, and 40.8 MPa (molality based) for CO 2 , H 2 S, and CH 4 , respectively. The Henry constants for these gases in ChClEg are calculated as 3.8 MPa, 1.4 MPa, and 60.1 MPa, using the charge scaling factor of 0.8. For ChClEg with a charge scaling of 0.8, the Henry constant of CO 2 is in agreement with the experimental value by Mirza et al. [63], but it is smaller than the value by Leron et al. [94] (Table 1). Therefore, the improvement in the accuracy of the computed Henry constant using the lower charge scaling factor depends on the considered experimental data (similarly for ChClU). The charge scaling value of 0.7 for ChClU results in an enhanced agreement of the simulation results with the experimental data of Xie et al. [79], and Liu et al. [60] (Table 1). Nevertheless, the differences between the computed Henry constants and the experimental data are still significant for ChClU. As discussed,  lower charge scaling factors result in lower densities and may therefore compromise the accuracy of the computed densities. Here, the reduced charge scaling factors resulted in densities of 1160 kg m −3 and 1070 kg m −3 for ChClU and ChClEg, respectively, which are lower than the computed values with the larger scaling factors (Fig. 1). The relative differences between the computed densities using different ionic charge scaling factors are therefore 2.1% and 2.8% for ChClU and ChClEg, respectively. Moreover, the lower charge scaling factors lead to an underestimation of experimental densities by 1.7% and 2.5%. Nevertheless, such density differences between simulation results and experimental data are not significant, and are also observed for the MD simulation results of Perkins et al. [37], and Doherty and Acevedo [39], using the GAFF and OPLS force fields, respectively, for ChClU. It should further be studied how these reduced charge scaling factors affect the liquid structure and other thermodynamic/transport properties (e.g., viscosity and diffusion coefficients) of the DESs. A closer agreement of the computed Henry constants of the gases in ChClU with experimental data of Xie et al. [79], and Liu et al. [60], however, will require even smaller ionic charge scaling factors. This will in turn result in less accurate calculations with respect to density and possibly other DES properties. It has been suggested that ionic charge scaling models can poorly predict the phase behavior of IL mixtures [124,127]. Cui et al. [128] investigated the effect of ionic charge scaling on solvation properties of ILs with respect to various solutes and concluded that such models under-predict the solute-solvent interactions and therefore free solvation energies, compared to experimental data. However, the authors stated that such differences are more pronounced in the case of more polar solutes (e.g., ammonia), whereas solute-solvent interactions are well reproduced by such models for non-polar solutes (e.g., CO 2 ) with dominant dispersion interactions. The authors further suggested that for the former case, explicit polarization of the models is necessary for accurate description of the interactions. It is also possible that DES properties can be more accurately computed in mixtures by using a polarizable force field, rather than the simple charge scaling method. It has been shown by García et al. [129] that the charge derivation scheme has a substantial influence on the computed structural and thermodynamic properties of choline chloride levulinic acid DES. Additionally, charge assignment can be carried out for either isolated ions/molecules or DES clusters. García et al. [129] argued that decisions regarding the charge assignment method should be made with great caution to result in accurate modeling of DESs. Therefore, other charge derivation schemes, for instance based on mixture properties, could improve the accuracy of solubility calculations. Recently, Schauperl et al. [130] developed a new version of the restrained electrostatic potential (RESP) method [131], RESP2, where partial atomic charges are derived based on a combination of gas-and aqueous-phase charges. In contrast, the original RESP model [131], by which Perkins et al. [37,38] computed the partial atomic charges of DESs (the charges used in this work) only performs gasphase quantum mechanical calculations. It is therefore interesting to see if new force field models developed based on mixture parameters, such as RESP2 [130], could lead to more accurate solubility computations. Additionally, the effect of LJ parameters can be examined on the computed properties of DES mixtures. For instance, one may consider using different mixing rules or modifying the used mixing rules to adjust the strength of interactions between specific unlike atoms [132][133][134]. In this work, the force field parameters were used as proposed in the original publications, and no further re-parameterization was performed for obtaining a closer agreement of the simulation results with experimental data. Although such re-parameterization may result in more accurate values for some of the Henry constants, there is no guarantee that the obtained parameters are transferable to other DES mixtures.

Conclusions
The Henry constants of various gases in ChClEg and ChClU DESs were computed using CFCMC simulations in the NPT ensemble. GAFF force field parameters were used for both DESs, in addition to OPLS parameters for ChClEg. The results were corroborated by additional isotherm calculations of CO 2 absorption in ChClEg in the osmotic ensemble. Densities and RDFs of the pure DESs were computed using MC simulations without applying the CFCMC technique. It was observed that the MC simulations could accurately reproduce the densities and RDFs of the DESs. The solubilities of the gases were, however, computed significantly lower (higher Henry constants) than experimental data reported in literature. The OPLS force field resulted in larger Henry constants of the gases in ChClEg (except for H 2 S) compared to the GAFF force field, suggesting that the GAFF force field may be more suitable for the solubility computations. The order of solubilities in both ChClU and ChClEg (using both force fields) at 328 K was H 2 S > CO 2 > CH 4 > H 2 > CO > N 2 , which is in reasonable agreement with experimental measurements available in literature. It was shown that smaller ionic charge scaling factors increase the solubilities of the gases in the DESs, but decrease the computed densities. Nevertheless, the scarcity and contradiction of experimental data raises uncertainty over the precise solubility values and the accuracy of predictions from the simulations. It is possible that the force field parameters used for the DESs in this study are not suitable for the solubility calculations, in combination with the force field parameters of the solute molecules. It is therefore proposed that force fields with explicit polarization terms or parameterized using mixtures of DESs with various compounds are examined for the computation of solubilities in DESs.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.