Choline-Based Ionic Liquids-Incorporated IRMOF-1 for H 2 S / CH 4 Capture: Insight from Molecular Dynamics Simulation

: The removal of H 2 S and CH 4 from natural gas is crucial as H 2 S causes environmental contamination, corrodes the gas stream pipelines, and decreases the feedstock for industrial productions. Many scientiﬁc researches have shown that the metal-organic framework (MOF) / ionic liquids (ILs) have great potential as alternative adsorbents to capture H 2 S. In this work, molecular dynamics (MD) simulation was carried out to determine the stability of ILs / IRMOF-1 as well as to study the solubility of H 2 S and CH 4 gases in this ILs / IRMOF-1 hybrid material. Three choline-based ILs were incorporated into IRMOF-1 with di ﬀ erent ratios of 0.4, 0.8, and 1.2% w / w , respectively, in which the most stable choline-based ILs / IRMOF-1 composite was analysed for H 2 S / CH 4 solubility selectivity. Among the three choline-based ILs / IRMOF-1, [Chl] [SCN] / IRMOF-1 shows the most stable incorporation. However, the increment of ILs loaded in the IRMOF-1 signiﬁcantly reduced the stability of the hybrid due to the crowding e ﬀ ect. Solvation free energy was then computed to determine the solubility of H 2 S and CH 4 in the [Chl] [SCN] / IRMOF-1. H 2 S showed higher solubility compared to CH 4 , where its solubility declined with the increase of choline-based IL loading. writing—review and editing, K.J., M.F.T., M.D.H.W., and M.A.; visualization, K.J., M.F.T., M.A.I.I., and M.N.N.; supervision, K.J., M.F.T., and M.D.H.W.; funding acquisition, K.J. and M.A.


Introduction
Hydrogen sulphide (H 2 S) is a highly toxic, colourless, flammable gas with a pungent odour. It is a gas that is often produced in the processing of crude oil and natural gases as it occurs in the subsurface as free gas and in oil due to its high solubility. H 2 S is also a by-product of the breakdown of organic waste as well as other industrial processes. It is significant to perform the removal of H 2 S as it causes operational problems in both oil and gas industries. It is highly corrosive to pipeline throughout the production facilities [1], leading to an economical deleterious effect because a washing plant is required to remove H 2 S in order to prevent corrosion and render the residual gas safe for domestic combustion. H 2 S also has adverse effects on the environment and human health. Exposure to high concentrations of H 2 S is extremely hazardous and can be fatal. Thus, the reduction of H 2 S emission becomes one of the main interests [2], especially in the biogas, syngas, and natural gas purification processes [3,4].
For the past decades, amine gas has been commonly used as a treating agent to separate H 2 S from natural gas in oil and gas industries. Single and binary mixtures of alkanolamines, such as diethanolamine (DEA), monoethanolamine (MEA), N-methyldiethanolamine (MDEA), and di-isopropanolamine (DIPA), are commonly used as the treating agents to remove H 2 S and CO 2 from natural gas [5]. However, it is necessary to find an alternative technique that is cleaner, with low energy consumption and cost, and at the same time is able to maintain or increase the production rate. The use of single and binary mixtures of alkanolamines as treating agents for the removal of H 2 S have some drawbacks in terms of performance and efficiency. Among them are the degradation of alkanolamine, the production of a corrosive side product due to oxidation reaction, and loss of alkanolamine inside the gas stream [6]. Plus, energy demand for the process is huge due to the wide energy gap on the reaction between the acid gases and alkanolamines, as well as the use of water as a solvent, which exhibits high heat capacity [7].
Ionic liquids (ILs) are currently used as a new technology to separate H 2 S from natural gases due to their high thermal stability, low volatility [8], and alteration of their functionalized group [9,10]. ILs are made up of cation and anion that exist in the form of molten salt or liquid with a melting point below 100 • C [11] or close to the ambient condition [12]. Due to these fascinating properties, ILs are widely used as solvents in numerous applications, including gas adsorption [13,14], catalysis and extraction [12,15], batteries and electrolyte [16][17][18], and molecular sensors [19]. Previous studies reveal the efficiency and capability of ILs to absorb H 2 S with numerous modifications and improvement. As the solubility of H 2 S in ILs is highly dependent on the partial pressure of H 2 S, most ILs can be classified as "physical" solvents and operated efficiently in the presence of H 2 S under high partial pressure and high concentration [20]. It has been reported that ILs with moderate basicity demonstrated better selectivity because the strong alkali groups interact excellently with H 2 S, making ILs a promising aspirant to separate H 2 S from other gases, such CO 2 and CH 4 .
Huang et al. [21] studied the performance of H 2 S absorption by using 1-alkyl-3-methylimidazolium carboxylate ILs at different temperatures, where the absorption capacities of H 2 S intensified compared to other common ILs. This is due to the alkalinity of the anions and also a small contribution from the alkyl chain lengths of the cations. Meanwhile, 1-hexyl-3-methylimidazolium chloride has been effectively used for the capture and conversion of H 2 S by using the Claus reaction, where H 2 S was converted to S 8 in a high conversion ratio by more than 96% in just 3 min [10]. Imidazolium-based ILs are widely used as liquid solvent for numerous applications, since they demonstrate excellent performance in gas adsorption and separation. However, the toxicity of imidazolium-based ILs has opened a new path towards finding a safe and environmentally-friendly IL. One of these paths are choline-based ILs, which have high a potential for practical use, since they are highly biodegradable and cost-effective. However, the efficiency of choline-based ILs in gas separation are hindered due to their high viscosity. Supporting ILs (SILs) are one of the promising solutions that are able to reduce the viscosity of ILs, thus improving the adsorption efficiency. The introduction of a metal-organic framework (MOF) as a support material will help in the immobilization of ILs, as they are dispersed all over the surface while preserving the unique properties of the MOF itself. In this concept, immobilization of ILs inside the MOF will reduce the strong association of ions, which contributes to their high viscosity, where the ions will be tethered on the solid surface. Therefore, SILs will be beneficial to both the MOF and ILs by improving the functionality and selectivity of the MOF towards gases, reducing the usage of ILs depending on the loading ratio, as well as increasing the gas solubility due to the presence of active sites.
MOF appears as one of the best materials to be used as the supporting absorber due to its intriguing features, such as high surface contact, high porosity [22][23][24][25][26] and pore volume, large gas storage, and high tunability [27][28][29], making it a suitable material for gas separation and adsorption. It is widely used in gas separation and storage [30][31][32], ion exchange [33,34], catalysis, and conductivity [35,36]. The MOF has been proposed as a SIL to improve the separation efficiency and enhance the physical nature of ILs by reducing the viscosity and increasing the gas solubility. The integration of both the MOF and ILs into a hybrid will reinforce the adsorption efficiency and, consequently, raise the H 2 S uptake. The proportion of ILs should be taken into consideration, since the adsorption capacity is affected by this factor. For example, the adsorption of a sulphurous compound (BT) was unsatisfactory with 33% loading of BMIM-Cl into MIL-101 [37]. This was due to the fulfilment of pore sites with BMIM-Cl in the MOF that hinder the adsorption availability. A similar observation was reported by Qiao et al. [38], where the adsorption of CO 2 was reduced by the addition of the amine group inside the MOF due to the unavailability of free pore volume. Selection of the MOF should be focussed on its ability to capture a high volume of guest molecules and on the selectivity of the desired guest molecules. Zheng et al. [39] stated that the design and synthesis of flexible MOFs, which can interact with certain guest molecules in a switchable way, would be appealing to explore. Furthermore, Abroshan and Kim [40] studied the stability of IRMOF-1 and IRMOF-10 after the incorporation of imidazolium and pyridinium-based ILs by molecular dynamics. They found that both IRMOF-1s are stable at low ILs loading but eventually collapse and deform at high IL loading, specifically at 15 ion pairs. Moreover, the stability of the IRMOFs is governed by the selection of different ILs, since each type contributes various stabilities towards the framework. Furthermore, Gupta et al. [41] conducted a simulation study on the adsorption of CO 2 /N 2 gas mixture inside the IL/IRMOF-1, where the CO 2 showed dominant interaction towards the anions. Meanwhile, Li et al. [42] computationally studied the solubility of H 2 S and CH 4 in IL/MOF composite using four types of imidazolium-based ILs, where Cu-TPDPAT was used as the MOF. It showed that the adsorption of H 2 S and its selectivity over CH 4 was improved upon IL loading compared to pristine Cu-TDPAT, in which the anion interacts dominantly with the gas compared to cation.
At present, computational study is utilised as a complement to experimental findings because it provides a detailed interaction and dynamic properties of the proposed study. Therefore, molecular dynamics (MD) simulation was conducted to provide an insight on the interaction, where IRMOF-1 was chosen as a potential MOF due to its high surface area and pore volume compared to other kinds of MOFs, specifically the zeolitic imidazolate framework (ZIF) and Cu-BTC. These remarkable properties help in providing a wide medium and capacity for gas accommodation. Meanwhile, choline-based ILs were selected as potential solvents due to their biodegradable, biocompatibility, and cost-effective qualities, compared to other common ILs. The combination of IRMOF-1 and choline-based ILs as a hybrid material will not only improve the gas adsorption performance but also promote the environmentally-friendly aspect, compared to common imidazolium and pyridinium-based IL/MOF. This work focussed on the prediction towards the stability of the IL/IRMOF-1 composite upon ILs incorporation, as well as the solubility of H 2 S and CH 4 inside the hybrid composite.  Figure 1) were packed inside the IRMOF-1 structure by using PyMOL. IRMOF-1 was obtained from the ChemTube 3D, originating from University of Liverpool [43]. The IRMOF-1 framework was built up by coordination of four oxide zinc atoms bonded to six carboxylate benzene rings in the subsequent arrangement to form a uniform porous pore. The optimized potential of liquid simulations (OPLS) force field was used, and the force field parameters were obtained from Chen et al., [44]. A Berendsen thermostat and Berendsen barostat were used for temperature and pressure control, respectively, in which the temperature was set up at 298 K, while the pressure was at 1 atm. A 2.6 × 2.6 × 2.6 nm cubic box was used to create the initial system of choline-based ILs/IRMOF-1. The IRMOF-1 was placed at the centre of the box and incorporated with three different weight percentages (W ILs/IRMOF1 0.4, 0.8, and 1.2%). Table 1 shows the number of ILs with respect to its ratios. For example, 14, 27, and 41 ion pairs of [Chl][SCN] were infused inside the IRMOF-1, where each number of ions pairs corresponded to a 0.4, 0.8, and 1.2% IL ratio. In this work, only one subunit of IRMOF-1 was used. The choline-based ILs and IRMOF-1 were packed together by using Packmol. Energy minimization was the first step in the simulation flow. In this step, there were two main algorithms; the steepest descent and the conjugate gradient scheme. A total of 100,000 steps were applied with a time step of 0.2 fs for each minimization step. This step was carried out to optimize the orientation of IL molecules in the system. The pre-production NVT run was applied for proper pre-equilibrium as final checking before the real simulation. The production simulation was performed in a conical ensemble (NPT) for 50 ns.

Simulation
Processes 2020, 8, x FOR PEER REVIEW 4 of 13 each number of ions pairs corresponded to a 0.4, 0.8, and 1.2% IL ratio. In this work, only one subunit of IRMOF-1 was used. The choline-based ILs and IRMOF-1 were packed together by using Packmol. Energy minimization was the first step in the simulation flow. In this step, there were two main algorithms; the steepest descent and the conjugate gradient scheme. A total of 100,000 steps were applied with a time step of 0.2 fs for each minimization step. This step was carried out to optimize the orientation of IL molecules in the system. The pre-production NVT run was applied for proper pre-equilibrium as final checking before the real simulation. The production simulation was performed in a conical ensemble (NPT) for 50 ns.

Solvation Free Energy
Solubility of H 2 S and CH 4 was predicted through solvation free energy simulation and calculation was done by using the Bennet acceptance ratio (BAR) method. Solvation free energy was computed with the same parameters and conditions as in the choline-based ILs/IRMOF-1 with an additional free energy parameter. A total of 21 sets of λ were used, starting from λ = 0 to λ = 1, with a total of 5 ns of simulation time for each λ. The energy of the systems was minimized using the steepest descent, followed by the conjugate gradient, with a total of 5000 steps for each algorithm. Then, the system was continued at 298 K and 1.0 bar for 100 ps in NVT, followed by 5 ns in NPT. The free energy difference between two states, i and j, was calculated using Equation (1), adopted from Bennett [45], where kB and T are the Boltzmann constant and temperature, respectively. Both Hi and Hj denoted the Hamiltonians at states i and j, while C was determined in repetition in order to satisfy the f(Hi − Hj + C)j = f(Hi − Hj − C)i condition. The f function was denoted from the Fermi function, as shown in Equation (2).

Solvation Free Energy
Solubility of H 2 S and CH 4 was predicted through solvation free energy simulation and calculation was done by using the Bennet acceptance ratio (BAR) method. Solvation free energy was computed with the same parameters and conditions as in the choline-based ILs/IRMOF-1 with an additional free energy parameter. A total of 21 sets of λ were used, starting from λ = 0 to λ = 1, with a total of 5 ns of simulation time for each λ. The energy of the systems was minimized using the steepest descent, followed by the conjugate gradient, with a total of 5000 steps for each algorithm. Then, the system was continued at 298 K and 1.0 bar for 100 ps in NVT, followed by 5 ns in NPT. The free energy difference between two states, i and j, was calculated using Equation (1), adopted from Bennett [45], where k B and T are the Boltzmann constant and temperature, respectively. Both H i and H j denoted the Hamiltonians at states i and j, while C was determined in repetition in order to satisfy the f ( condition. The f function was denoted from the Fermi function, as shown in Equation (2).
The free energy difference was calculated from Equations (3) and (4), where N j and N i indicate the number of coordinate frames at λ i and λ j , respectively. Technically, the forward (FW) and backward (BW) energy differences should be overlapped adequately to achieve the convergence of this process, which is the crucial aspect for the BAR method [45]. The data would then be reliable if the overlap integral relied minimally on 0.01, as proposed by Bruckner and Boresch [46], in which the overlapping could be observed from the normalized histograms of the energy differences. Otherwise, the repetition of the process would not be converged or unreliable data would be generated due to the insufficient overlap.
In Equations (5) and (6), the λ symbol represents the λ-value once the trajectory is generated. The overlap could be detected by overlaying the histograms from both Equations (5) and (6) with respect to the ∆U l FW and ∆U l+1 BW values. However, the actual amount of overlap cannot be obtained directly through the visual inspection. Hence, the overlap integral OI l,l+1 BAR was introduced as the numerical guide to provide the proper point for improvement in the existing free energy protocol. The interval λ l , λ l+1 is formulated as Equation (7), where h b (X) represents the normalized histogram of quantity X with a total number of B bins. The non-zero value of OI l,l+1 BAR was contributed from the configuration in the overlap zone, in which 0 indicates that there is no overlap zone, while 1 represents a perfect overlap.
The solubility of each gas is calculated using Equation (8) [47].

Stability of ILs/IRMOF-1
After the incorporation of choline-based ILs into IRMOF-1, the root mean square deviation (RMSD) was calculated to determine the stability and sustainability of choline-based IL/IRMOF-1 composites. Table 2 shows the RMSD value for ILs/IRMOF-1 at three different weight ratios. IRMOF-1 is well known for its fragility and poor stability in the aqueous phase, where it tends to collapse upon contact with water [48]. With the presence of water molecules in pristine IRMOF-1, the calculated RMSD was 1.423 nm at 298 K. This high deviation justified the instability of IRMOF-1, but upon incorporation of ILs inside the IRMOF-1, its stability significantly improved by preventing direct contact with water. Therefore, it can be deduced that the incorporation of ILs inside the IRMOF-1 caused some modification in molecular structure from the original state. Based on Table 2 Moreover, the stability of the hybrid material varied, according to the number of ILs inside the IRMOF-1. Based on Table 2, the calculated RMSD value increased as the weight percentage of ILs multiplied. This result indicates that a small number of ILs was sufficient to maintain the stability of the composite. This phenomenon could be explained as follows: By increasing the number of ILs inside the MOF, the excess number of ILs created a crowding effect which hindered the interionic interaction of the attractive ions towards the IRMOF-1. Consequently, the interaction between choline-based ILs and the IRMOF-1 deteriorated, thus reducing the stability of the choline-based ILs/IRMOF-1.

Interaction of Choline-Based ILsİnside the IRMOF-1
From the RMSD value, the infusion of choline-based ILs inside the IRMOF-1 affected the stability of the framework, in which the interaction of the cation and anion of the ILs also varied upon their incorporation. In order to quantify the interaction between the anion and cation inside the IRMOF-1 framework, the radial distribution function (RDF) was plotted, as shown in Figure 2. The interaction of the opposite ions observed was lower in the bulk phase, as shown in Figure 2a, compared to the hybrid materials phase. Meanwhile, their interaction significantly increased inside the IRMOF-1 framework, in which the highest interaction between opposite ions was exhibited by 0.4% w/w ([Chl][SCN]/IRMOF-1) loading. The strong interaction could be explained by the confinement effect being posed by the IRMOF-1 framework, in which the ions inside the framework were packed closely together. In addition, the different charges posed by [Chl] + and [SCN] − promoted strong electrostatic forces, as can be seen by the red peak in Figure 2a.
of the opposite ions observed was lower in the bulk phase, as shown in Figure 2a, compared to the hybrid materials phase. Meanwhile, their interaction significantly increased inside the IRMOF-1 framework, in which the highest interaction between opposite ions was exhibited by 0.4% w/w ([Chl][SCN]/IRMOF-1) loading. The strong interaction could be explained by the confinement effect being posed by the IRMOF-1 framework, in which the ions inside the framework were packed closely together. In addition, the different charges posed by [Chl] + and [SCN] − promoted strong electrostatic forces, as can be seen by the red peak in Figure 2a.   packed inside the IRMOF-1, which shows the crowding effect by [Chl][SCN] at a high percentage of ILs (refer to Figure 3b,c). The accumulation of [Chl] + and [SCN] − ions inside the IRMOF-1 framework with various IL loadings was inspected in 2D visualisation by mapping through the density map taken from the z-axis direction, as shown in Figure 4. It can be observed that the pore volume of IRMOF-1 was fully occupied by [Chl] + as it was larger in size compared to [SCN] − (Figure 4b). As shown by Figure 4c, a low density of [SCN] − accumulated inside the pore due to its small size.

Solubility and Solubility Selectivity of H 2 S/CH 4 Inside [Chl][SCN]/IRMOF-1
The solubility of H 2 S and CH 4 inside the [Chl][SCN]/IRMOF-1 was calculated from the solvation free energy. Table 3 shows the excess chemical potential (µ ex ) and calculated solubility of each gas inside the hybrid material for three different loadings of [Chl] [SCN]. The solubility rate of the gases was quantified by the value of µ ex , in which the solubility of each gas intensified as the negative value of µ ex increased. Based on Table 3, H2S and CH4 in the pure ILs and IRMOF-1 exhibited lower solubility compared to those in IL/IRMOF-1 composites, respectively. The µ ex value of H2S inside the pristine IRMOF-1 was 0.122 kcal/mol, suggesting its poor solubility. This indicates that IRMOF-1 had a weak interaction with H2S, which reflects its low affinity and selectivity in capturing gas. As for the pure ILs, we predict that the viscosity was the limitation for the gas solubility. This is because viscous IL possesses strong association between its anion and cation, which will hinder interaction with gases. molecules qwas less favourable compared to H2S, due to the non-polarizability of CH 4 [50].

Solubility and Solubility Selectivity of H 2 S/CH 4 Inside [Chl][SCN]/IRMOF-1
The solubility of H 2 S and CH 4 inside the [Chl][SCN]/IRMOF-1 was calculated from the solvation free energy. Table 3 shows the excess chemical potential (µ ex ) and calculated solubility of each gas inside the hybrid material for three different loadings of [Chl] [SCN]. The solubility rate of the gases was quantified by the value of µ ex , in which the solubility of each gas intensified as the negative value of µ ex increased. Based on Table 3, H 2 S and CH 4 in the pure ILs and IRMOF-1 exhibited lower solubility compared to those in IL/IRMOF-1 composites, respectively. The µ ex value of H 2 S inside the pristine IRMOF-1 was 0.122 kcal/mol, suggesting its poor solubility. This indicates that IRMOF-1 had a weak interaction with H 2 S, which reflects its low affinity and selectivity in capturing gas. As for the pure ILs, we predict that the viscosity was the limitation for the gas solubility. This is because viscous IL possesses strong association between its anion and cation, which will hinder interaction with gases. Meanwhile, H 2 S and CH 4 showed a negative value of µ ex in IL/IRMOF-1, which indicates that both gases were soluble in the composite. However, H 2 S dominated in the solubility compared to CH 4 [49]. On the other hand, the interaction of this IL with CH 4 molecules qwas less favourable compared to H 2 S, due to the non-polarizability of CH 4 [50].  [51].
In order to determine the binding and interaction of each gas toward [Chl] + and [SCN] − inside the IRMOF-1, the RDF was plotted. Figure [52]. The highest solubility of H 2 S was observed at 303 K with the µ ex value of −1.116 kcal mol −1 and 0.266 mol L −1 bar −1 at 14 IL (0.4% w/w) loading.  Furthermore, the calculated solubility of H 2 S at 14 IL (0.4% w/w) loading decreased from 0.266 to 0.184 mol L −1 bar −1 at 298 K and 333 K, respectively. This trend was similar to 27 IL (0.8% w/w) and 41 IL (1.2% w/w) loading for both H 2 S and CH 4 . Since anions contributed a major effect towards the solubility of the gas [42,53], the rise in temperature significantly affected the interaction between anions and gases. The decline in gas solubility as the temperature rose could be explained from the weakening of the bonding in the intermolecular interaction between IL and gas. This is because [SCN]existed as a conjugate base and acted as a hydrogen-bond acceptor [54]. Therefore, [SCN] − was able to form hydrogen bonding with the gases, which in turn allowed the gases, especially H 2 S, to solubilise [55]. However, the increase in temperature weakened the bonding and reduced the intermolecular interaction of IL-gas. This is in agreement with other experimental works [56,57].

Conclusions
In conclusion, the stability of the choline-based ILs/IRMOF-1 composite was successfully predicted by MD simulation. [