Accelerate the Electrolyte Perturbed-Chain Statistical Associating Fluid Theory–Density Functional Theory Calculation With the Chebyshev Pseudo-Spectral Collocation Method. Part II. Spherical Geometry and Anderson Mixing

To improve the efficiency of electrolyte perturbed-chain statistical associating fluid theory–density functional theory (ePC-SAFT-DFT) calculation of the confined system, in this work, first, the Chebyshev pseudo-spectral collocation method was extended to the spherical pores. Second, it was combined with the Anderson mixing algorithm to accelerate the iterative process. The results show that the Anderson mixing algorithm can reduce the computation time significantly. Finally, based on the accelerated ePC-SAFT-DFT program, a systematic study of the effects of the temperature, pressure, pore size, and pore shape on the CO2 solubilities in the ionic liquids (ILs) confined inside the silica nanopores was conducted. Based on the simulation results, to obtain high CO2 solubilities in the ILs confined in silica, a better option is to use the silica material with a narrow spherical pore, and the IL-anion should be selected specifically considering that it has a more significant impact on the absorption enhancement effect.


INTRODUCTION
Mitigating CO 2 emission from fossil-fueled power plants as well as from transports has become an urgent and worldwide research topic, in which CO 2 separation is often needed (MacDowell et al., 2010;Boot-Handford et al., 2014). Ionic liquids (ILs) are promising absorbents for CO 2 separation due to their extremely low vapor pressure, high CO 2 solubility, as well as low-energy usage for solvent regeneration (Brennecke and Gurkan, 2010;Ramdin et al., 2012;Zhang et al., 2012). However, the viscosity of pure ILs is relatively high compared with common organic solvents, causing a significant decrease in the mass and heat-transfer rates. Using supported ILs has been proposed as a promising solution for practical applications. This can take the advantage of high gas selectivity in ILs, and also, the high surface area of the supported materials can reduce the impact of high viscosity, improve the gas transfer, and hence increase the absorption rate (Zhang et al., 2006;Ren et al., 2012;Romanos et al., 2014).
Research has been conducted to address the confinement effect on the gas solubility in ILs via experiments and molecular simulations (Baltus et al., 2005;Ilconich et al., 2007;Ilconich et al., 2007;Zhang et al., 2010;Iarikov et al., 2011;Ren et al., 2012;Banu et al., 2013;Shi and Luebke, 2013;Romanos et al., 2014;Santos et al., 2014;Budhathoki et al., 2017;Shen and Hung, 2017). According to previous research, several factors will affect the CO 2 solubility inside the confined ILs, for example, temperature, pressure, pore size, and shape of porous materials (Baltus et al., 2005;Zhang et al., 2006;Ilconich et al., 2007;Zhang et al., 2010;Iarikov et al., 2011;Ren et al., 2012;Banu et al., 2013;Shi and Luebke, 2013;Romanos et al., 2014;Santos et al., 2014;Budhathoki et al., 2017;Shen and Hung, 2017). However, to screen a suitable IL, optimizing the structure of supported material and operation conditions by experiment or molecular simulations is time and cost consuming, considering the fact that the huge number (10 18 ) of possible ILs can be synthesized (Wasserscheid and Thomas, 2008), as well as the wide temperature and/or pressure range in applications. Therefore, it is desirable to develop a theoretical model to predict the properties of confined IL-CO 2 systems.
The classical density functional theory (DFT) is considered as an efficient theoretical method for studying the confined properties (Tripathi and Chapman, 2005;Qiao et al., 2018;Shen et al., 2013;Xu et al., 2008;Sauer and Gross, 2017;Camacho Vergara et al., 2019). In addition, in our previous work (Ji et al., 2012;Ji et al., 2014;Shen et al., 2015;Ji and Held, 2016;Sun et al., 2019), electrolyte-perturbed-chain statistical associating fluid theory (ePC-SAFT) (Cameretti et al., 2005) has been developed to represent the thermodynamic properties of IL systems. Moreover, the developed ePC-SAFT has been combined with DFT (ePC-SAFT-DFT) to describe the properties of IL and CO 2 /IL confined in nanopores with acceptable results (Shen et al., 2018). Recently, in order to calculate the properties of the confined ILs with ePC-SAFT-DFT efficiently, the Chebyshev pseudo-spectral collocation method (Yatsyshin et al., 2012;Nold et al., 2017) was implemented to accelerate the ePC-SAFT-DFT calculation . However, only the slit-shaped and cylindrical pores have been considered previously, while for the spherical cavity, the corresponding method has not been available. In addition, an advanced iteration method is required for replacing the simple Picard iteration to accelerate ePC-SAFT-DFT calculation further. Anderson mixing is an elaborate iteration method that has been used in the work by Mairhofer and Gross (2017) and Shen et al. (2021), showing desirable performance in accelerating DFT computing. Therefore, replacing the Picard iteration with Anderson mixing can be an effective strategy.
In this work, the Chebyshev pseudo-spectral collocation method was extended to the spherical geometry, where an expression of 9-3 Lennard-Jones potential for the spherical cavity was derived. In addition, Anderson mixing was used to accelerate the ePC-SAFT-DFT calculation further. Based on the modified program, the CO 2 solubility of ILs confined in the silica nanopores was chosen as the representative to conduct the investigation, considering silica is a promising supporting material for ILs (Shi and Luebke, 2013), and the effects of temperature, pressure, pore structures, and IL-ions were investigated systematically.

Electrolyte Perturbed-Chain Statistical Associating Fluid Theory-Density Functional Theory
According to DFT, in the presence of a solid surface, the grand potential Ω at grand canonical ensemble is given by the equation: where A is the Helmholtz free energy, ρ i (r) is the molecular density of component i at position r, μ i is the chemical potential, m i is the number of segments in a chain for component i, and V i,ext (r′) is the nonelectrostatic external field acting on the segment of component i. The Helmholtz free energy A in ePC-SAFT-DFT for a system can be expressed as (Qiao et al., 2018;Shen et al., 2018): where A id is the ideal free energy, A hs , A chain , A disp , and A ion are the excess free energies due to hard-sphere repulsions, chain connectivity, dispersive, and electrostatic interactions, respectively. The model performance has been verified in a previous work (Shen et al., 2018), where the model predictions in the density profiles of ionic fluids in the charged pores were compared with the molecular simulation results with high consistency. The details of A hs , A chain , and A disp have been described elsewhere (A hs : Barker and Henderson, 1967;Rosenfeld, 1989;Roth et al., 2002;Yu and Wu, 2002a;Yu and Wu, 2002b;A chain : Tripathi and Chapman, 2005; A disp : Shen et al., 2013). The A ion is composed of two terms, the Coulomb term (A col ) and the Debye-Hu€ckel term (A DH ). A DH accounts for the short-range electrostatic interaction caused by the inhomogeneous distribution at the nearby region of ions, which is based on the Debye-Hu€ckel (DH) theory. A col accounts for the Helmholtz free energy caused by the long-range electrostatic interaction due to the unsymmetrical distribution of cation and anion. The expressions of A col and A DH can be found in Li and Wu (2006) and Shen et al. (2018), respectively. All these terms are involved in the ePC-SAFT-DFT calculations for the IL-CO 2 systems.
Minimization of the grand potential with respect to the density profile of component i yields the following Euler-Lagrange equation: where q i is the charge of component i, and ψ(r) is the mean electric potential at position r, which can be obtained by solving the Poisson's equation (Li and Wu, 2006). The expressions of functional derivatives δA c [ρ(r)] δρ i (r) are described in Rosenfeld et al. (1997), Gross (2009), Roth (2010, and Sun et al. (2021).

Evaluation of the Convolution-Like Integrals for Spherical Cavity
The expressions of the functional derivatives δA c [ρ(r)] δρ i (r) are described in Rosenfeld et al. (1997), Gross (2009), Roth (2010, and Sun et al. (2021). All the convolution-like integrals involved in ePC-SAFT-DFT calculation can be classified as: where R c represents the weighting distances, which are: In Eq. 4b, θ is the Heaviside step function, and δ is the Dirac delta function. The σ i and d i in Eq. 4b represent the temperature-independent and -dependent hard sphere segment diameter (Barker and Henderson, 1967), respectively. The d i,DH is the weighting distance for the Debye-Hu€ckel term, which can be referred to in Shen et al. (2018). The λ was set as 1.5 in order to be consistent with the bulk PC-SAFT model (Gross and Sadowski, 2001). The convolution integrals in Eq. 7a map f(r′) to I(r) with these integral kernels.
For spherical cavities, f(r′) and I(r) only vary with the radial direction, and the analytic expressions of I(r) reads (Groh and Mulder, 2000):

Evaluation of Mean Electric Potential Distribution Inside Spherical Cavity
In the spherical symmetric distribution, the Poisson's equation reduces to: where q(r) is the total charge at position r.
For a spherical cavity with diameter R, by solving Eq. 6 with the boundary conditions ψ(R) ψ w and dψ(0) dr 0, the expression for the mean electric potential reads: According to Eq. 7, the boundary electric potential ψ w is required, which corresponds to a specific charge density on the internal surface of a nanopore. This value can be obtained via "trial and error" until the charge of ions in the nanopore is equal in magnitude to that with the opposite sign on the internal surface of the nanopore.ü

Chebyshev Pseudo-Spectral Collocation Method
The Chebyshev pseudo-spectral collocation method for DFT modeling was developed by Yatsyshin et al. (2012). In the Chebyshev pseudo-spectral collocation method, for onedimensional DFT calculation with an N-point discretization scheme, the density or the weighted density profile over the whole computation domain is determined from the density or the weighted density at a prescribed set of collocation points {z k }, k = 1,2 . . . , N using the barycentric form (Baltensperger, 2002;Berrut and Trefethen, 2004): where the primes indicate that the first and last terms in the sums are divided by 2. The collocation points (z k ) can be obtained by a conformal map from the Chebyshev collocation points (x k ) (Baltensperger, 2002;Berrut and Trefethen, 2004): The integrals associated with the DFT calculation can be evaluated with the Clenshaw-Curtis quadrature (Clenshaw and Curtis, 1960).

Implementation of Chebyshev Pseudo-Spectral Collocation Method in Electrolyte Perturbed-Chain Statistical Associating Fluid Theory-Density Functional Theory for Spherical Cavity
In the ePC-SAFT-DFT calculation, three domains need to be discretized for interpolation. The first domain is for the density profile, the second one is for the weighted density function profiles in the hard-sphere term, and the third one is for the weighted density function profile in the Debye-Hu€ckel term.
For the spherical cavity with diameter R, we defined a diameter R′ as: The density profile is considered in the domain (0, R') based on the coordinate system illustrated in Figure 1. In other words, this domain is discretized for interpolating the density profile function. In general, the density profile vibrates dramatically near the wall, which implies that more collocation points are required in these regions compared with the middle of the nanopore. In this work, the conformal map proposed by Bayliss and Turkel was used to map the Chebyshev collocation points {x k } to domain (0, R') (Bayliss and Turkel, 1992): The λ and μ in Eq. 11a read: where: Equation 11a maps the Chebyshev points in the transformed coordinate into the points that cluster in the physical coordinate near R′ 2 (β + 1) with a density that is large when α is large. In this work, the β was set as: where d is the mean diameter of all components. Therefore, the collocation points are clustered in the physical coordinate near (R′ − d 2 ), where high peaks have occurred in the nearby region. The α was set as 2 in this work. Similar approaches were used for the discretization of another two domains, i.e., the FIGURE 1 | Schematic diagram of the domains considered in a spherical cavity.

Evaluation of Convolution-Like Integral
As pointed out by Yatsyshin et al. (2012), the maps (Eq. 4) involved in DFT calculation can be carried out by matrix-vector products with discrete data based on interpolation (Eq. 8) and Clenshaw-Curtis quadrature. In the spherical cavities, for a map with original space discretized into N 2 points and the image space discretized into N 1 points, the map can be represented by N 1 × N 2 matrix (Yatsyshin et al., 2012): where a i and b i are the lower and upper limits of integral for the ith discrete point, W is the corresponding integral kernel in Eq. 5, and y qi is the qth Chebyshev grid of the ith discrete point in the image space: w q is the corresponding Clenshaw-Curtis weight, and the values at M Chebyshev grids are used to evaluate each integral. In this work, M was set as 45.
For the singularities in the expression of Eq. 13b (i.e., y qi z k ), the corresponding elements can be evaluated based on the following expression :

Evaluation of the Mean Electric Potential
Equation 7 can be represented with matrix-vector products with the discrete data: where ψ and q are composed of ψ and q at the position of collocation points, respectively. The two N × N matrixes D 1 and D 2 read: where L′ and U′ are composed of: where L and U are the lower and upper triangular matrices with all elements being unity, respectively. The elements in the (N-1) × N matrix Q 1 read: The elements in (N-1) × N Q 2 read: y qi can be evaluated with Eq. 13c with the a i and b i defined in Eq. 14h. The singularities in equations Eq. 14e and Eq. 14g can also be evaluated based on Eq. 13d. Q 1 and Q 2 are the linear operators implementing the piecewise integrations between two adjacent collocation points, while L′ and U′ are used to sum up the results of these piecewise integrations for obtaining the final integration results of Eq. 7.

Solving Euler-Lagrange Equation With Anderson Mixing
In the Chebyshev pseudo-spectral collocation method, the Euler-Lagrange equation (Eq. 3) is transformed into a vector equation that can be solved numerically. In our previous work, the Picard iteration was used. In this work, the Anderson mixing is implemented to solve Eq. 3 iteratively (Anderson, 1965). In the Anderson mixing, the Euler-Lagrange equation can be rewritten as: where k represents the Boltzmann constant, ρ i,b represents the bulk phase density of component i, and the superscript res represents the residual quantities.
Eq. 15a can be solved iteratively by (Anderson, 1965): where S is the relaxing factor, m k min(m, k), m was set as 50 in this work, the subscript i represents the ith component, and γ k j is determined by: The constrained optimization can be transformed to unconstrained optimization (Fang and Saad, 2008;Walker and Ni, 2011): The successive optimization (Eq. 15d) can be solved efficiently by the updating QR factorization, and the necessary Matlab code is given in Walker (2011). In order to avoid divergence, Ns steps of Picard iterations can be performed at the beginning and then switched to the Anderson mixing procedure. SCHEME 1 | Calculating the density profile of ionic liquid (IL)-CO 2 system inside nanopores with specific surface charge Q w using the general scheme combined with the Anderson mixing.

General Scheme Combined With Anderson Mixing
A general scheme for calculating the density profile of confined systems with ionic contribution was proposed in Sun et al. (2021), as illustrated in Scheme 1. In a loop of Scheme 1, the Euler-Lagrange equation needs to be solved twice. For the first time, solving the Euler-Lagrange equation in a loop, Ns (steps of Picard iteration performed before the Anderson mixing procedure) can be set between 1,200 and 1,500, while for the second time, due to the initial guess of density profile being not far away from the equilibrium density profile, we do not perform Picard iterations prior to the Anderson mixing procedure. (i.e., Ns = 0).
In Scheme 1, T, p, μ i , and x i refer to those in the bulk phase, which were obtained with ePC-SAFT. H (or R) and Q w are the features of nanopores. ρ in represents the initial guess of the density profile, and ρ e represents the calculated equilibrium density profile with the boundary electric potential ψ w . The charge absorbed per area of nanopore Q can be computed with: where the amount of component i adsorbed per surface area ρ i for the slit-like pores can be evaluated by: for the cylindrical pores: and for the spherical cavities:

Model Ionic Liquid-CO 2 Confined in Nanopores
Following ePC-SAFT for the ILs in the bulk systems (Ji et al., 2012), an ionic liquid molecule is composed of one IL cation and one IL anion. Each individual IL ion was modeled as a nonspherical species with repulsion, dispersive attraction, and Coulomb interactions. The 9-3 Lennard-Jones potential was used to represent the nonelectrostatic interaction between silica and fluid. The nanopore was modeled as an infinitely large slit, infinitely long cylinder, or spherical cavity. The 9-3 Lennard-Jones potential for large slit and infinitely long cylinder have already been presented in the literature (Fitzgerald et al., 2003;Siderius and Gelb, 2011;Lee, 2016), while for a spherical cavity with diameter R, the 9-3 Lennard-Jones potential can be obtained from: where r is the distance of the fluid molecule from the center of the spherical cavity, ρ atom is the solid atom density, σ si is the  effective solid-fluid diameter, and ε si is the potential representing the interaction between surface and fluid segment. σ si and ε si can be determined with the Berthelot-Lorentz combining rules: where σ s and ε s are the size and potential parameters of a solid surface, respectively, and ε i is the potential parameters of fluid segment i. Here, the used potential parameters for silica are σ s 3.0 Å, ε s 0.8 kJ/mol, and 4πρ atom 0.5/Å 3 (Pinilla et al., 2005). The ePC-SAFT parameters for ILs were taken from the previous work (Ji et al., 2014), while those for CO 2 were taken from Gross and Sadowski (2001). For the electroneutral silica nanopore, the amount of cation and anion adsorbed per surface area should be equal. The solubility of CO 2 in the confined IL with a neutral surface is defined by: where ρ IL ρ IL−cation ρ IL−anion , according to the charge neutrality condition.

Efficiency of the General Scheme Combined With Anderson Mixing
Calculating the density profile of [C 6 mim] [Tf 2 N]-CO 2 confined in electronic neutral silica pore with different structures at 333 K and 16.1 bar was selected as an example here to demonstrate the performance of the general scheme combined with the Anderson mixing, and the calculation efficiency was compared with the general scheme only using the Picard iterations. In the calculation, the width of the slitshaped pore and the diameters for the cylindrical pore and spherical cavity are 5 nm. We used 120 collocation points to represent the density profile inside the slit-shaped pore (functional derivatives only need to be calculated in half of the collocation points due to the symmetry of the density profile), while 60 collocation points were used for the cylindrical pore and spherical cavity. The relaxing parameters used in all these three cases are 0.001.
The calculations were performed on a computer with an AMD core Ryzen 7 PRO 4750U and X64 processor. The version of the compiler is Matlab 2018b. Figure 2 compares the time required for calculations.
As illustrated in Figure 2, using the Anderson mixing can improve the calculation efficiency significantly. The time needed for calculating the matrices required in the Chebyshev pseudo-spectral collocation method is listed in Table 1.
As listed in Table 1, the time required for calculating the matrices for the slit-shaped pore and spherical cavity can be ignored. However, the time for the cylindrical pore is pronounced. The reason is that too much time is used for evaluating the Legendre complete elliptic integral of the first and second kinds. This may be a problem in massive ePC-SAFT-DFT calculations, for example, adjusting the model parameters from experimental data, in which the equilibrium density profile needs to be solved multiple times.
However, for systems at the same temperature and pore diameter (cylindrical pore), the elements of these matrices are the same . Therefore, for the sake of saving computation time, when modeling systems confined in the cylindrical pore, these matrices can be used repeatedly for the systems at the same temperature and pore diameter.

Model CO 2 Solubilities of Confined Ionic Liquids
Based on the efficient algorithm discussed above, ePC-SAFT-DFT can be used in a wide range of IL systems to obtain results efficiently. In this work, several [C n mim]-based IL-CO 2 systems confined in silica nanopore were selected to investigate the effects of temperature, pressure, IL ions, as well as the size and shape of the pore. As the nanopores of different silica materials have been roughly assumed as slit-like (Li et al., 2005;Yang and Yue, 2007), cylindrical Beck et al., 1992;Maddox et al., 1997), and spherical (Nandiyanto et al., 2009) in the previous theoretical work, in this work, these three pore-shape models were adopted.

General View of CO 2 Confined In Silica Nanopores
The calculated density profile of [C 6 mim] [Tf 2 N]-CO 2 confined in 25 Å slit silica pore at 10 bar and 323.15 K is presented in Figure 3.
As pointed out by Ho et al. (2013), two competitive mechanisms affect the solubility in a nanopore: One is the adsorption near the pore surface, and the other is the absorption inside the IL. These can be seen in the density profile of CO 2 . According to the results illustrated in Figure 3, the peak in the density profile near the surface of nanopores corresponds to the first mechanism, while in the middle of the nanopore, the density profile of CO 2 tends to the bulk density, which is the reflection of the second mechanism. The absorption enhancement in the confined ILs is mainly contributed by the first mechanism. In this work, the additional solubility x a is used to identify the absorption enhancement: x a x confined − x bulk (19)

The Influence of Temperature
In general, with the temperature increase, the CO 2 solubility in the bulk ILs will decrease. According to our calculation results, CO 2 solubility in [C 6 mim] [Tf 2 N] confined in slit-shaped pore also decreases with the increase in temperature. This is consistent with the observation by Mirzaei et al. (2017). Typical examples are shown in Figures 3A and 4A. However, with the increase in temperature, the solubility of CO 2 in [C 6 mim] [Tf 2 N] confined in SiO 2 decreases greater than that of the bulk IL. For example, at 1 bar, when the temperature increase from 298.15 K to 373.15 K, the solubility of CO 2 in the bulk [C 6 mim] [Tf 2 N] decreases by  Frontiers in Chemistry | www.frontiersin.org January 2022 | Volume 9 | Article 801551 11 about 0.02 mol CO 2 /mol IL. Under the same condition, the calculated solubilities of CO 2 in the confined [C 6 mim] [Tf 2 N] reduce by about 0.03 mol CO 2 /mol IL. It indicates that confined ILs may make it easier to desorb CO 2 .
In order to study the effect of temperature on the absorption enhancement effect, we analyzed the relation between the additional solubility and temperatures. As illustrated in Figures 4B and 5B, the additional solubility decreases with increasing temperature. In order to interpret this, the calculated density profiles of [C 6 mim][Tf 2 N]-CO 2 confined in 25 Å slit silica pore at 10 bar and different temperatures are illustrated in Figure 6. To quantitatively describe the effect of the first mechanism, the ratios of the CO 2 density at the first peak near the pore surface (which is also the maximum density inside the nanopore) to its bulk density in ILs were calculated as illustrated in Figure 7. According to the results shown in Figure 7, the ratio decreases with the increase in the temperature when the pressure keeps constant, and consequently, the additional solubility decreases. The same phenomenon can be observed at other pressures and in other pore structures, as listed in Supplementary Table S1.

The Influence of Pressure
Most commercial ILs are physical absorbents for CO 2 . In general, the CO 2 solubility in ILs will increase significantly with increasing pressure. Figure 8A illustrates the CO 2 solubility at 298.15 K and different pressures in the 25 Å slit-shaped pore. From Figure 8B, The CO 2 solubilities in confined ILs also increase with the increase in the pressure. In addition, the increase in CO 2 solubility in the confined [C 6 mim] [Tf 2 N] is more significant than that of the bulk [C 6 mim] [Tf 2 N] under the same condition. For instance, when the pressure increases from 1 to 50 bar, the solubility of CO 2 in the confined [C 6 mim] [Tf 2 N] increases by about 0.78 mol CO 2 /mol IL, while that in bulk [C 6 mim] [Tf 2 N] increases by about 0.70 mol CO 2 /mol IL.
The additional solubility also increases with increasing pressure, as illustrated in Figure 8B. To have a deep insight into this observation, the density profiles of [C 6 mim] [Tf 2 N]-CO 2 confined in 25 Å slit silica pore at 298.15 K and different pressures illustrated in Figure 9 were taken as one example for the detailed analysis. According to the results shown in Figure 9, with the increase in pressure, the adsorption of CO 2 near the pore surface increases, while the adsorption of IL ions near the pore surface decreases. This can be interpreted that at high pressure, a large amount of CO 2 inside the nanopore leads to a stronger competitive adsorption capacity of CO 2 in the pore surface than that of IL ions. Consequently, the additional solubility increases with the increase in pressure. The same phenomenon can be observed under other conditions, as listed in Supplementary Table S1.

The Influence of Pore Size and Shape
The CO 2 solubilities of confined ILs will also be affected by the pore size and shape. According to the calculation results, the adsorption will be enhanced more significantly in the smaller nanopore, which is consistent with the results of Shi and Luebke (2013). Two typical examples are shown in Figure 10. As illustrated in Figure 10, the solubilities of CO 2 increase with the decrease in the slit-shaped pore.
The effect of the shape of nanopores was also investigated. According to our calculation results, for nanopores with the same size (i.e., the width for slit-shaped pore and the diameter for cylindrical pore and spherical cavity), ILs confined in the spherical cavity have the highest CO 2 solubilities, while those confined in the slit-shaped pore have the lowest. Typical examples are illustrated in Figure 11. This indicates that using supported material with spherical nanopores may lead to better absorption capacity. The calculated density profiles of the three shapes at 298.15 K and 30 bar are presented in Figure 12. It shows that the pore shape affects the density profile near the surface of the pore significantly. For IL ions, the density at the first peak follows the trend that ρ (slit-shaped) > ρ (cylindrical) > ρ (spherical), while for CO 2 , the opposite trend can be observed. This is because the inward curved surfaces impede the accumulation of large molecules (e.g., IL ions) near the surface. Therefore, under the same situation, the CO 2 solubilities (x) follow the trend that x (spherical) > x (cylindrical) > x (slit-shaped).

The Influence of Cation
In most cases, the IL-cations are chain-like substances. Therefore, in this part, the effect of alkyl-chain length was studied. Following this, ePC-SAFT-DFT was used to model the CO 2 solubilities of [C 4 mim][Tf 2 N], [C 6 mim] [Tf 2 N], and [C 8 mim] [Tf 2 N] inside silica nanopores at different temperatures, pressures, pore widths, and pore shapes. The calculated results are listed in Supplementary Table S1.
As illustrated in Figures 13A and 14A, for the ILs in the same homologous series, the CO 2 solubilities in the confined ILs increase with increasing alkyl-chain length in cation, which is  Frontiers in Chemistry | www.frontiersin.org January 2022 | Volume 9 | Article 801551 13 the same as that in the bulk ILs. According to Figure 13B and Figure 14B, at low pressures, the additional solubility also increases with the increase in the alkyl-chain length. However, this is not true at high pressures, as illustrated in Figure 14B. This indicates that the absorption enhancement of ILs with longer alkyl-chain decreases faster with the increasing pressure than the ILs with shorter alkyl-chain.
As demonstrated in Figure 15, with the increase in pore width, the additional solubility decreases, especially for the IL with longer alkyl-chain length in the IL cation. The additional solubilities of confined ILs in different pore shapes are illustrated in Figure 16. In the spherical cavity, the additional solubilities increase more significantly with increasing alkyl-chain length in IL cation than that in the slit-shaped and cylindrical pores.

The Influence of Anion
In the bulk phase, the IL anion affects the CO 2 solubilities more significantly than IL cation (Anthony et al., 2005). According to   In Figure 18, the calculated additional solubilities of several confined ILs in different pore shapes are presented. The absorption enhancement of anion is more significant in the spherical cavity as shown in Figure 18. For example, in the slitshaped pore at 30 bar, the deviation of the calculated additional solubilities between [C 4  ], the additional solubilities are considerably higher than those of [C 6 mim] [Tf 2 N] in the spherical cavity at high pressures. This indicates that changing the IL anion of ILs confined in the spherical cavity may be promising to obtain a large absorption enhancement.
Up to here of this section, the effects of temperature, pressure, IL ions, as well as the size and shape of pores on confined CO 2 solubilities were investigated. In order to further improve model performance toward practical applications, the ePC-SAFT-DFT model will be combined with different pore models (pore shape and pore size distribution) with the parameters adjusted from experimental data in our future work based on the newly measured experimental data.

CONCLUSION
In this work, the Chebyshev pseudo-spectral collocation method was combined with the Anderson mixing algorithm to further accelerate the ePC-SAFT-DFT calculation. The results show that the computing time can be significantly reduced with the Anderson mixing algorithm. This makes the ePC-SAFT-DFT model more effective in screening promising ILs. However, calculating matrices used in the Chebyshev pseudo-spectral collocation method for the cylindrical pores requires a certain computing time, while in a massive computation, these matrices can be reused to save computation time.
Using the ePC-SAF-DFT model, the CO 2 solubilities of ionic liquid confined in silica nanopores were studied. The results show that the CO 2 solubilities of confined ILs are always higher than that in bulk ILs under the same condition, and the absorption enhancement effect will be significantly affected by pressure, pore widths, pore shapes, and IL anion. Based on the simulation results, to obtain high CO 2 solubilities in silica-confined ILs, a better option is to use silica   Frontiers in Chemistry | www.frontiersin.org January 2022 | Volume 9 | Article 801551 material with a narrow spherical pore. In addition, the IL anion should be selected specifically considering that the category of IL anion has a significant impact on the absorption enhancement effect.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding authors.