Molecular screening of ionic liquids for CO₂ absorption and molecular dynamic simulation

2.1 Screened ILs and COSMO-RS calculations

ILs can theoretically be designed by combining different anions and cations; 25 cations and 20 anions are selected to combine 500 IL solvents. Their structures and names are listed in Tables S1 and S2 of the supporting information. The cations include pyrazolium, phosphonium, ammonium, imidazolium, pyridinium, pyrrolidinium, guanidinium, and piperidinium. The anions include acetate, phosphate, sulfate, and dicyanamide. COSMOtherm was used to predict the solubility of CO₂ in the above-mentioned ILs. The standard operating procedure of COSMOtherm contains two steps. First, quantum calculation of the involved molecule species is performed by Gaussian or Turbmole suit program. Second, statistical thermodynamics calculation is carried out within COSMOtherm. All screened IL structures could be seen as the combination of ion pairs, namely, cation and anion. The structures of ion pairs, which had been early optimized using four different theoretical levels and incorporated into the program, were directly loaded from the database. As suggested in the study by ref. [25], a molecular model of counter ions was applied in COSMOtherm calculations. ILs are treated as equimolar mixtures of cations and anions. The parameterization (BP_TZVPD_FINE_20.ctd) contains the intrinsic parameters for the calculation and is used. The σ-profile, which contains the screening charges on the molecular surface, was used as an input in the COSMOtherm to calculate the solubility of CO₂ in ILs at near ambient pressure and temperature.

2.2 Molecular dynamics simulations

MD simulations were carried out using GROMACS 2018.2 program package [26]. The force field, General AMBER Force Field (GAFF), was used in the MD simulations. Each anion or cation in ILs was geometrically optimized at B3LYP/6-311 + G(d,p) in Gaussian 16 package [27]. The atom RESP charges of ILs were calculated using Multiwfn 3.7 program [28,29]. The bonded and non-bonded parameters used in the simulation are obtained by ACPYPE online server [30]. A leap-frog algorithm was used to integrate Newton’s equations of motion, and the time step of 2fs was adopted in all the simulations. The temperature was controlled by a velocity rescaling thermostat with a coupling constant of 0.4 ps [31]. The pressure was maintained by Parrinello–Rahman barostat with a coupling constant of 2.0 ps [32]. The periodic boundary conditions were applied in all directions. The bonds were constrained with the LINCS algorithm [33]. The cut-off for the short-range electrostatics and Lennard–Jones interactions were both set to 1.0 nm, and the particle mesh Ewald algorithm [34] was used to deal with the long-range electrostatic interactions. For each system, energy minimization of the simulation system was first performed with the steepest descent algorithm to generate an appropriate starting configuration and then equilibrated for 1 ns at isothermal and fixed volume (NVT) ensemble to reach the reference temperature of 298.15 K and the reference pressure of 1 bar, respectively. Then, 20 ns trajectory was produced at the isothermal and isobaric (NPT) ensemble, and the last 10 ns trajectory was used for statistical analysis.

Free energy of CO₂ in pure water and ILs was calculated using the BAR (Bennett Acceptance Ratio) method with 21 equidistant sets of λ values from λ = 0 to λ = 1. TIP4P water model was used. Coordinates, velocities, and δH/δλ were saved to disk at every ten steps for the post-processing of BAR. Each node was subjected to energy minimization with 20,000 steps of the steepest descent method. The simulation was then pre-equilibrated at NVT ensemble with no position restrains. After NVT equilibration, the process was continued with 2 ns simulations at NPT at 298 K for all systems. Data collection was done from the last 1 ns.

3.1 CO₂ solubility in ILs

The solubility of CO₂ in 500 ILs is calculated at 298 K and 100 kPa using the method mentioned above, and the results are tabulated in Table S3. Figure 1 shows the trend and difference of CO₂ solubility in the screened ILs.

Figure 1

CO₂ relative solubility in ILs calculated by COSMOtherm.

The previous studies have investigated cation and anion effects on CO₂ solubility [35,36,37]. The combination of cation and anion can affect CO₂ solubility. Generally, CO₂ solubility in ILs strongly depends on the choice of the anion. Aki et al. also studied the influence of the anion with the same 1-butyl-3-methylimidazolium cation [11]. They proved that the fluorination of the anion can improve CO₂ solubility. The trend of CO₂ solubility in ILs predicted by COSMOtherm is consistent with the previous work. It was concluded from Figure 1 that ILs with B(Hex)₃P cation have much higher CO₂ solubility than any other ILs, and CO₂ solubility in ILs with bFAP anion was also much higher than any other ILs. CO₂ solubility in various ILs is much different. Among the 500 screened ILs, the highest CO₂ solubility is 0.084 mole fraction in [B(Hex)₃P][bFAP], while the lowest CO₂ solubility is 0.002 mole fraction in [MPrAm][NO₃] at 298 K and 100 kPa.

The previous experiments explored ILs containing nitrile anion that had a higher CO₂ solubility [38,39]. Their work showed that CO₂ solubility in [emim][B(CN)₄] (1-ethyl-3-methyl imidazolium tetracyanoborate) is 0.027 mole fraction. Compared with that, ILs containing [bFAP] anion have much higher CO₂ solubility.

3.2 MD simulations of ILs and ILs/CO₂

A few of ILs were selected according to the CO₂ solubility data. Specifically, bFAP anion was paired with four different cations, namely, [EMPyr](1-ethyl-2-methylpyrazolium), [B(Hex)₃P](butyl-trihexyl-phosphonium), [(Me)₅isobuGua](n,n,n,n,n-pentamethyl-n-isopropyl-guanidinium), and [BEIM](1-butyl-3-ethyl-imidazolium). Two kinds of MD simulations were performed. The first set was to examine CO₂ solvation free energy in ILs. The second was to mimic the equilibrium and dynamic properties of pure ILs and the IL/CO₂ system. For the calculation of CO₂ solvation free energy in IL, 1 CO₂ molecule was put into the center of the cubic box. For the pure IL simulation system, a cubic box of 8 nm × 8 nm × 8 nm was created with 256 pairs of cation and anion. For IL/CO₂ system, the simulation box was constructed by 10 CO₂ molecules and 256 ion pairs. The initial configurations of the MD simulations were built by the Packmol package [40].

4.1 Henry constant and solvation free energy of CO₂ in ILs

Henry constant, as another measure of CO₂ solubility in ILs, could also be computed using COSMOtherm. The computed Henry constants of CO₂ in pure water and the four selected ILs are listed in Table 1. Additionally, to understand and compare the CO₂ absorption ability in ILs and water from the point of solvation free energy, the CO₂ solvation free energy in pure water was also computed by the BAR method and COSMOtherm. The values of the solvation free energy of CO₂ (abbreviated as ΔG) in pure water and four different ILs are listed in Table 1. The calculated results using the BAR method are 3.80 kJ mol⁻¹ in pure water, −10.18 kJ mol⁻¹ in [EMPyr][bFAP], −11.37 kJ mol⁻¹ in [BEIM][bFAP], −12.21 kJ mol⁻¹ in [(Me)₅isobuGua][bFAP], and −10.43 kJ mol⁻¹ in [B(Hex)₃P][bFAP], respectively. Compared with the solvation free energy of CO₂ in pure water, the negative values clearly showed that the CO₂ molecule is more easily dissolved in IL than in water.

Though the four ILs could dissolve CO₂ more easily than water, it was concluded that there were different CO₂ solubility capabilities among the four selected ILs. Though the values of ΔG obtained by COSMOtherm were much smaller than that obtained by the BAR method, they all showed the same trend. The solvation free energy of CO₂ in studied ILs followed the increasing order: [EMPyr][bFAP] < [BEIM][bFAP] < [(Me)₅isobuGua][bFAP] < [B(Hex)₃P][bFAP], which agrees with the relative CO₂ solubility predicted by COSMOtherm.

4.2 Structural properties

Some important structural properties of the simulated systems were investigated through MD simulations. The density, solvent-accessible area (SAS), and free volume (FV) of all simulated systems are summarized in Table 2.

The maximum density corresponds to [EMPyr][bFAP] (1.844 g cm⁻³), and the minimum density corresponds to [B(Hex)₃P][bFAP] (1.431 g cm⁻³). It was noticed that the density of the IL/CO₂ binary system was a little smaller than that of the corresponding pure IL. FV within ILs was correlated to the gas solubility. CO₂ molecules were preferentially solvated within the vicinity of the anions [41]. Figure 2 shows the representative vicinity within ILs. It was noticed that the four ILs with bFAP anion all had a larger FV of more than 35%. The larger the FV was, the more CO₂ molecules were solvated. It was consistence with the CO₂ solubility in the four ILs.

Figure 2

Representative snapshot of the FV space (green) within bFAP-based ILs: (a) [BEIM][bFAP], (b) [B(Hex)₃P][bFAP], (c) [EMPyr][bFAP], and (d)[(Me)₅isobuGua][bFAP].

4.3 Radial distribution function of ILs and ILs/CO₂ system

To study the influence of CO₂ absorption on the overall IL structure, the equilibrium structures of ILs and CO₂ are presented by the radial distribution function g(r). The g(r) of the constituent compounds are calculated to study cation–anion, cation–cation, and anion–anion interactions in IL and shown from a to c in Figure 3. The g(r) of cation, anion, and CO₂ in the ILs/CO₂ system is also calculated and shown from (d) to (i) in Figure 3.

Figure 3

Radial distribution functions in ILs and IL/CO₂ systems at 298 K and 100 kPa: (a) cation–anion in IL, (b) cation–cation in IL, (c) anion–anion in IL, (d) cation–anion in IL/CO₂, (e) cation–cation in IL/CO₂, (f) anion–anion in IL/CO₂, (g) CO₂–cation in IL/CO₂, (h) CO₂–anion in IL/CO₂, and (i) CO₂–CO₂ in IL/CO₂.

The g(r) provided deep insights into the interactions that occurred in the studied system. It was revealed from Figure 3(a–c) that cation–anion interaction of both ILs is stronger than those of cation–cation or anion–anion interactions. Figure 3a shows the interaction between the cation and anion in ILs. Similar observations were found that the weaker interaction between the cation and the anion affected the CO₂ solubility in ILs [42]. Taking [B(hex)₃][bFAP] as an example, the intense cation–anion g(r) peak appears at about 0.67 nm. The g(r) peak of cation–anion in [BEIM][bFAP] was about 0.56 nm, while in Figure 3b and c, the cation–cation and anion–anion g(r) of [B(hex)₃][bFAP] peaks appear at 1.20 and 1.10 nm, respectively. It was clear in Figure 3d that the cation–anion g(r) shape in IL/CO₂ systems did not change obviously compared with pure IL systems. It was consistent with the previous simulation studies which showed that the addition of CO₂ molecules does not change the structures of ILs [36,43]. The interactions between CO₂ and ILs ions are shown in Figure 3(g and h). The CO₂–anion interactions in the four IL/systems all display g(r) peaks at 0.54 nm. It is obvious that CO₂ molecules are interacting stronger with the anions as compared with the cations since the CO₂–anion peak is located at a lower distance than that of the CO₂–cation peak.

The electrostatic and the van der Waals energies between ions and CO₂ interactions were also computed from the MD simulation systems and are listed in Table 3.

For IL and IL/CO₂ systems, the total interaction energies of cation–cation and anion–anion mainly showed unfavorable interactions. As for cation–anion, it was favorable interaction and electrostatic interaction made more contributions to the total interaction than Van der Waals interaction. Electrostatic interaction of cation and CO₂ was larger than that of anion and CO₂. Van der Waal’s interaction of anion and CO₂ was larger than that of cation and CO₂. It showed that the interaction of cation and CO₂ was mainly electrostatic interaction, while the interaction of anion and CO₂ was mainly Van der Waals interaction.

4.4 Mean-squared displacement (MSD) and Diffusion coefficients

The dynamic properties of ILs and IL/CO₂ systems are quantified using MSD and diffusivity. Figure 4(a–e) show the MSDs of cation, anion, and CO₂ in ILs and ILs/CO₂ at 298 K.

Figure 4

Comparison of the mean square displacement of of cations, anions, and CO₂ in IL and IL/CO₂ at 298 K: (a) [BEIM] and [bFAP] in pure [BEIM][bFAP], (b) [B(Hex)₃P] and [bFAP] in pure [B(Hex)₃P][bFAP], (c) [EMPyr] and [bFAP] in pure [EMPyr][bFAP], (d) [(Me)₅isobuGua] and [bFAP] in pure [(Me)₅isobuGua][bFAP], and (e) CO₂ in the four ILs.

The initial rapid increase in MSD within 1 ns showed that a longer simulation time was required to evaluate the self-diffusion coefficients of ions or CO₂ in both IL and IL/CO₂ systems. The slopes of the MSDs for the cations are steeper than those for the anions in IL and IL/CO₂ systems. It was possibly due to the fact that a larger anion size resulted in a slower movement than the cation.

The diffusion of cation, anion, and CO₂ in ILs and ILs/CO₂ systems has also been determined. The diffusion coefficients of the constituent compounds were obtained by the linear fitting of the slope of MSD. Table 4 reports the diffusivities of cations, anions, and CO₂ in the system.

The diffusion of cation, anion, and CO₂ in ILs and ILs/CO₂ systems have also been determined. The diffusion coefficients of the constituent compounds were obtained by the linear fitting of the slope of MSD. Table 4 reports the diffusivities of cations, anions, and CO₂ in the system. When CO₂ molecules were dissolved in the IL, it was synergistic competition. It is clear from Table 4 that the diffusion coefficient of CO₂ was much higher than those of cations and anions. It was mainly due to the small size of CO₂ molecules. For CO₂ molecule in IL/CO₂ system, the self-diffusion coefficient followed the decreasing order: [B(Hex)₃P][bFAP]/CO₂ > [(Me)₅isobuGua][bFAP]/CO₂ > [BEIM][bFAP]/CO₂ > [EMPyr][bFAP]/CO₂, which agreed with the predicted CO₂ solubility in ILs.