Gas Adsorption and Diffusion Behaviors in Interfacial Systems Composed of a Polymer of Intrinsic Microporosity and Amorphous Silica: A Molecular Simulation Study

We investigate the adsorption and diffusion behaviors of CO2, CH4, and N2 in interfacial systems composed of a polymer of intrinsic microporosity (PIM-1) and amorphous silica using grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulations. We build model systems of mixed matrix membranes (MMMs) with PIM-1 chains sandwiched between silica surfaces. Gas adsorption analysis using GCMC simulations shows that gas molecules are preferentially adsorbed in microcavities distributed near silica surfaces, resulting in an increase in the solubility coefficients of CO2, CH4, and N2 compared to bulk PIM-1. In contrast, diffusion coefficients obtained from MD simulations and then calibrated using the dual-mode sorption model show different tendencies depending on gas species: CO2 diffusivity decreases in MMMs compared to PIM-1, whereas CH4 and N2 diffusivities increase. These differences are attributed to competing effects of silica surfaces: the emergence of larger pores as a result of chain packing disruption, which enhances gas diffusion, and a quadrupole-dipole interaction between gas molecules and silica surface hydroxyl groups, which retards gas diffusion. The former has a greater impact on CH4 and N2 diffusivities, whereas the latter has a greater impact on CO2 diffusivity due to the strong quadrupole-dipole interaction between CO2 and surface hydroxyls. These findings add to our understanding of gas adsorption and diffusion behaviors in the vicinity of PIM-1/silica interfaces, which are unobtainable in experimental studies.


INTRODUCTION
Carbon dioxide (CO2) capture and separation technologies are active research areas for addressing the rapid increase in atmospheric CO2 concentration causing severe environmental issues such as an acceleration of global warming, 1 sea-level rise, 2 and ocean acidification. 3 Compared to conventional yet prevalent amine-based methods, 4,5 membrane-based gas separation technologies 6,7 have several advantages: high-energy efficiency, small footprint, operational simplicity, and modularity. Specifically, separation membranes employing microporous polymers, [8][9][10][11][12] defined as those with pore sizes smaller than 2 nm, are promising candidates for CO2 separation technologies because their pore sizes are comparable to the molecular dimensions of CO2 and other smaller gases. Furthermore, microporous polymers are more advantageous than other inorganic materials in terms of low cost, film-forming ability, processability, chemical diversity, and chemical and thermal stability.
Gas transport in separation membranes follows the solution-diffusion mechanism, 13 where the permeability ( ) is defined as a product of the diffusion coefficient ( ) and the solubility coefficient ( ), and the permselectivity of gas over gas is given as ⁄ . The current research efforts have been significantly devoted to surpassing the Robeson's upper bound 14 which describes the trade-off between permeability and permselectivity. To this end, various microporous polymers have been developed, including hyper-crosslinked polymers, 8 conjugated microporous polymers, 9 thermally rearranged polymers, 10 and polymers of intrinsic microporosity (PIMs). 11,12 Specifically, PIMs are a class of ladder-like polymers with rigid and contorted backbone structures. Their rigid contortion sites hinder the efficient packing of molecular chains, resulting in large free volumes, high surface areas, and unique pore structures, leading to increased CO2 permeability and high permselectivity. Since the inception of archetypal PIM-1, 15,16 various PIMs and their derivatives have been developed and exhibited performances exceeding the upper bound. 11 Mixed matrix membranes (MMMs) 17,18 formed by incorporating nonporous or porous nanofillers into polymer matrices have also attracted considerable attention as promising candidates to further push the upper bound. Incorporated nanofillers, whether porous or not, alter local polymer-chain packing and mobility, and porous nanofillers such as zeolites, zeolitic imidazolate frameworks (ZIFs), and metal-organic frameworks (MOFs) offer additional internal permeation paths. In either case, the structural properties at the polymer/nanofiller interface significantly differ from those in the bulk polymeric phase; hence gas adsorption and diffusion behaviors at the interface play a critical role in determining the separation performance in MMMs. PIM-1-based MMMs have been developed with various nanofillers such as PIM-1/silica, 19,20 PIM-1/MOFs, 21 PIM-1/ZIF-8, 22 PIM-1/multi-walled carbon nanotubes, 23 and PIM-1/porous aromatic frameworks. 24 Although improvements in gas permeability have been reported for these PIM-1-based MMMs, the molecular-level mechanisms of gas adsorption and diffusion at the polymer/nanofiller interface are yet to be fully explored. Molecular simulations are attractive techniques to probe molecular transport in interfacial systems at the atomistic scale.
Some molecular simulation studies have investigated the structural properties of PIM-1/ZIF-8 and PIM-EA-TB/ZIF-8 interfaces, 25,26 gas transport through a PIM-1/ZIF-8 interface, 27 and structure and gas transport at a polyimide/zeolite interface. 28 However, to the best of our knowledge, no simulation study has been reported that investigates gas adsorption and diffusion behaviors at the PIM-1/silica interface. Although silica is a nonporous material, surface modification of silica nanoparticles has been reported to enhance their dispersibility in the polymer phase and increase the performance of PIM-1-based MMMs, surpassing the upper bound. 29 We expect a molecular-level understanding of gas adsorption and diffusion behaviors at the PIM-1/silica interface to be a valuable starting point toward a rational design of silicaincorporated MMMs.
In this study, we investigate the adsorption and diffusion behaviors of CO2, CH4, and N2 in interfacial systems composed of PIM-1 and amorphous silica using grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulations. We build two types of model MMMs comprising PIM-1 chains sandwiched by silica surfaces, referred to as MMM-1 and MMM-2, for systems with a single chain and two chains, respectively. Silica surfaces significantly affect PIM-1 chain packing, which is more pronounced for MMM-1 with a shorter distance between the silica surfaces. Gas adsorption analysis using GCMC simulations shows that gas molecules are preferentially adsorbed in microcavities distributed near silica surfaces, increasing the solubility coefficients of CO2, CH4, and N2 relative to bulk PIM-1. Meanwhile, diffusion coefficients obtained from MD simulations and then calibrated using the dual-mode sorption model 30,31 show different tendencies: CO2 diffusivity decreases in MMM-1 and MMM-2 compared to PIM-1, while CH4 and N2 diffusivities increase. These disparities are due to competing effects of silica surfaces: the emergence of larger pores resulting from chain packing disruption, which enhances gas diffusion, and a quadrupole-dipole interaction between gas molecules and silica surface hydroxyl groups, which retards gas diffusion. More specifically, the former has a greater impact on CH4 and N2 diffusivities, whereas the latter has a greater impact on CO2 diffusivity due to CO2's strong quadrupole moment. These findings provide intriguing insights into gas adsorption and diffusion behaviors in the vicinity of PIM-1/silica interfaces, which are inaccessible to experiments.

Construction of PIM-1.
The construction of PIM-1 (Figure 1a) is performed using an in silico polymerization scheme as implemented in the Polymatic code. 32 This scheme allows for constructing an amorphous polymer by simulated polymerization of the constituent monomers, and its application has been proven for a range of polymers, including PIMs. [33][34][35] In this study, PIM-1 molecules are represented using a united-atom (UA) approach to reduce computational costs. The GAFF 36 and TraPPE-UA 37,38 force fields are used for the bonded and nonbonded interactions, where the Lorentz-Berthelot combination rules are adopted for the cross Lennard-Jones (LJ) interactions.
The atomic charges are derived from quantum chemical calculations based on HF/6-31G(d) level of theory, followed by a restrained electrostatic potential (RESP) charge fitting procedure. 39 This charge scheme is the default approach in GAFF parametrization and has been shown to be suitable for a wide range of organic molecules. 36 More details on the force field parameters and atomic charges can be found elsewhere. 40 The timestep is set to be 1 fs, and the cutoff radius for the LJ and short-range electrostatic interactions is set to be 15 Å with the long-range electrostatic interactions evaluated using a particle-particle particle-mesh method 41 with a precision factor of 10 −4 . All MD simulations are performed using the LAMMPS software package, 42 and the visualizations of MD systems are performed with OVITO. 43 Initially, 200 PIM-1 monomers are randomly packed in a cubic simulation box at low density (~ 0.3 g/cm 3 ) under the periodic boundary conditions. Each monomer contains predefined linker atoms that can participate in bond formation if the bonding criteria for a cutoff radius (6 Å) and participating monomer orientations are met (planarity and directionality).
Energy minimization and short MD simulations are performed throughout the process to avoid high-energy conformation. This simulated polymerization is repeated until all linker atoms are used, yielding a single PIM-1 chain with a degree of polymerization (DP) of 200.
In the present study, we prepare five structures of PIM-1 membranes to account for statistical uncertainties. The bulk density ( ), which is directly obtained from MD simulation, is related to experimentally obtained skeletal density ( ) as 48 where denotes the mass of the polymer and denotes the pore volume. The bulk density is calculated to be We also construct five PIM-1 structures with DP = 200 that are made up of eight molecular chains. The structure factor of the 8-chain membranes is almost identical to that of the 1-chain membranes, as shown in Figure S1 in Supporting Information, indicating that there are no apparent size effects. As a result, the discussion of PIM-1 membranes that follows is limited to the results obtained from the 1-chain systems.   Amorphous silica systems are constructed from β-cristobalite structure using the melt-quench method. 55 The β-cristobalite structure consists of 8232 atoms in a cubic simulation box of 50.12 × 50.12 × 50.12 Å 3 under the periodic boundary conditions ( Figure S2 in Supporting Information). Here, the interatomic interactions are represented by the Morse-style potential 56,57 with the cutoff radius of 10 Å, with the force field parameters provided in Table S1 in Supporting Information. Initially, the β-cristobalite structure is melted at 7000 K for 1 ns in the NVT ensemble. Next, the system is quenched to 300 K at a cooling rate of 4 K/ps, and intermediate constant-temperature simulations are performed for 100 ps at T = 6000, 5000, …, 1000, and 300 K to alleviate any stresses built up in the system due to the rapid cooling. 58 The cooling rate of 4 K/ps has also been used for constructing amorphous silica at room temperature in previous studies, 55,59 leading to realistic representation of its structural properties. Finally, the system is relaxed in the NPT ensemble for 200 ps at T = 300 K and P = 1 atm, followed by equilibration in the NVT ensemble for 200 ps at T = 300 K.
The amorphous silica structure thus obtained ( Figure S2 in Supporting Information) exhibits a density of 2.23 g/cm 3 , which agrees with previous studies. 59 respectively, which are used as cutoff radii for defining coordination numbers ( -and -).
Coordination numbers for the amorphous silica structure are provided in Tables S2 and S3 in Supporting Information.
To create silica surfaces, a vacuum layer is formed along the z-axis with a length sufficient to prevent interactions between top and bottom surfaces across the periodic boundary ( Figure   1c). The surfaces are relaxed in the NVT ensemble for 10 ps at T = 300 K. Then, the dangling silicon ( -4) and oxygen ( -2) atoms are saturated with hydroxyl groups and hydrogens. The resultant area density of silanol groups is 3.08 /nm 2 (Figure 1c), which is within the range of literature values (2.6-4.6 /nm 2 ). 60 The silica surfaces contain 3% Q 2 (i.e., two silanol groups per superficial silicon), 53% Q 3 (i.e., one silanol group per superficial silicon), and 44% Q 4 (i.e., siloxane bridges without silanol groups) environments. This surface chemistry corresponds to neutral silanol-terminated surfaces at pH values between 2 and 4 which previously underwent thermal treatment in experiment. 61 In reality, surface chemistry varies depending on the type of substrate, surface ionization relating to pH and particle size, synthesis protocol, and prethermal treatment. More specific models reflecting these details can be considered using the INTERFACE force field, 62 which should also be investigated in future work.
To model hydroxylated silica surfaces, the force field is switched to the CHARMM water contact angle (CWCA) force field, 63 which has been designed to reproduce experimental water contact angles on silica. The force field parameters and atomic charges 55 66 Under these conditions, an NVT-MD simulation is performed at 300 K for 400 ps using a Langevin thermostat 67 with a damping factor of 100 fs to equilibrate the surface silanol groups.

Construction of PIM-1/Silica Hybrid Systems.
We build five MMMs from a single PIM-1 chain (DP = 200) and amorphous silica, and five MMMs from two PIM-1 chains and amorphous silica; these are referred to as MMM-1 and MMM-2, respectively ( Figure 1d). In this section, we describe how to build a PIM-1/silica hybrid system for the case of MMM-1.
The PIM-1 chain obtained in Section 2.1 is first unwrapped in all directions and placed in an orthorhombic cell with the dimensions shown in Table S6 in Supporting Information. The system is then compressed and relaxed using the 21-step equilibration scheme (Table 1), resulting in cell lengths in the x and y directions ( ) around 50 Å (Table S6 in Supporting Information), which are very close to the cell lengths in the x and y directions of the silica slab system ( 49.71 Å) built in Section 2.2. Under the constant-volume condition, the PIM-1 system is slightly elongated in the z-direction to exactly match ( ) with ( ). The PIM-1 chain is then unwrapped in the z-direction and inserted into the vacuum layer of the silica slab system, as illustrated in Figure S4 in Supporting Information.
Finally, as proposed by Semino and coworkers, 25 the constructed PIM-1/silica system is compressed and relaxed in the z-direction using the 21-step equilibration scheme (Table 1), yielding the equilibrated system shown in Figure 1d. To maintain the silica structure, the silica framework composed of silicon and bridging oxygen atoms is kept rigid, and compression is performed only in the z-direction.
In the case of MMM-2, two PIM-1 chains are first placed in an orthorhombic cell with the dimensions shown in Table S7 in Supporting Information. The following steps are the same as for MMM-1.

Evaluation of Gas Solubility.
The solubility coefficient is evaluated from the slope of an adsorption isotherm in the dilute limit as where is the concentration of adsorbed gases and is the pressure.
In general, sorption of gases onto glassy polymers can be described using the dual-mode sorption (DMS) model 68-70 as where is Henry's law constant, is the Langmuir monolayer sorption capacity, and is the Langmuir affinity constant. The DMS model is based on two types of sorption sites subject to Henry's law dissolution and Langmuir-type sorption and has been shown to accurately represent the gas sorption in PIM-1. 71,72 This study evaluates the adsorption of CO2, CH4, and N2 in PIM-1, MMM-1, and MMM-2 through grand canonical Monte Carlo (GCMC) simulations using the LAMMPS software package. 42 GCMC simulations perform exchanges of molecules with an imaginary gas reservoir at the specified temperature ( ) and chemical potential ( ). 73 The gas molecules are modeled with the TraPPE-UA force field, 37,38 where CO2 and N2 are treated as rigid molecules, and CH4 is represented by the united-atom model with a single interaction site. GCMC simulations for a variety of chemical potentials are performed to obtain adsorption isotherms for each gas species in each membrane. GCMC simulations are run at 300 K in the μVT ensemble for each chemical potential (or pressure), with the membrane framework held rigid during gas adsorption. After each GCMC step, a 100-ps MD simulation is performed at 300 K in the NVT ensemble to equilibrate the system. This GCMC-MD cycle is repeated until the number of gas molecules in the system converges sufficiently. Before the gas adsorption simulations in the membranes, the relationship between the chemical potential and pressure is obtained for each gas species through bulk gas-phase GCMC simulations to take account of the nonideality of the gases at high pressures ( Figure S5 in Supporting Information).
We note that the present study does not consider plasticization and swelling of the membranes upon gas adsorption. It has been reported that gas adsorption at high pressures induces the swelling of PIM-1, 74-76 increasing the free volume and chain segment mobility.
Although some studies have proposed the methods for simulating the swelling of polymeric membranes, 77-79 rigid nature of the silica framework in the present study allows for the volume dilation only in the interface-normal direction (Figure 1d), prohibiting a proper evaluation of the membrane swelling behavior. Rather, the present study mainly focuses on the PIM-1/silica interface that leads to the disruption of polymer-chain packing and generation of microcavities in the vicinity of the interface; these factors significantly affect adsorption and diffusion behaviors of the gases in the membranes.

Evaluation of Gas Diffusivity.
The self-diffusion coefficient can be evaluated from MD simulation using the Einstein relationship 52 as where is the number of gas molecules, is the dimensionality of the diffusing molecules under consideration, and is the center-of-mass position of molecule at time . The term in brackets is the so-called mean-squared displacement (MSD).
To determine the self-diffusion coefficients of CO2, CH4, and N2 in PIM-1, MMM-1, and MMM-2 at dilute conditions, 10 gas molecules are randomly inserted into the pore space within the polymer region for each gas-membrane combination. The gas molecules' trajectories are then sampled in the NVE ensemble, while the temperatures of the membrane frameworks are kept at 300 K using the Nosé-Hoover thermostat. 44,45 The three-dimensional MSD of the gas molecules ( 3) is evaluated in PIM-1, while the two-dimensional MSDs in the xy-plane ( 2) are calculated in MMM-1 and MMM-2 because the PIM-1 regions are bounded by the nonporous silica frameworks (Figure 1d). We note that the self-, collective, and transport diffusion coefficients coincide in the limit of zero loading of gas molecules. 80,81 In the present study, inserting 10 gas molecules into each membrane corresponds to the gas pressures of at most 0.02, 0.2, and 0.6 atm for CO2, CH4, and N2, respectively (see Figure 4). For these small pressure values, the obtained self-diffusion coefficients would be good proxies for the transport diffusion coefficients. To confirm convergence to the zero-loading diffusion coefficient, it is necessary to calculate self-and transport diffusion coefficients using equilibrium and nonequilibrium MD simulations with varied gas loadings 82 and extrapolate the diffusivityloading curves to the zero-loading limit, which is beyond the scope of this study.
Unfortunately, diffusion coefficients calculated using MD simulations cannot be directly where denotes the diffusion coefficient associated with the Langmuir-type sorption and is the gas pressure upstream of the polymer film. The functional form of is given by 31 , ,

Gas Adsorption Behavior.
Figure 4a depicts the adsorption isotherms of CO2, CH4, and N2 in PIM-1. All of the isotherms have convex pressure curves, and CO2 is the most condensable, with an initial loading at low pressures that is significantly higher than those of CH4 and N2. The CO2, CH4, and N2 solubility coefficients in PIM-1 are calculated using eq 2 and tabulated in Table 2 Table 3.  isotherms, as shown in Figures 4c and 4d. MMM-1 has the highest solubility coefficient for each gas species, followed by MMM-2 and PIM-1, as shown in Table 2. To explain the increases in solubility coefficients in MMM-1 and MMM-2, gas density distributions in representative MMM-1 are evaluated for a range of gas pressures, as shown in Figure 5, and those in representative MMM-2 are shown in Figure S6 in Supporting Information. As shown in Figure   5a, for CO2 pressure < 1 atm, more CO2 molecules are adsorbed in the PIM-1/silica mixing regions (defined in Figure 2a   Finally, as shown in Figure 4, the adsorption isotherms of CO2, CH4, and N2 in PIM-1, MMM-1, and MMM-2 are fitted by the DMS model (eq 3) with the fitting parameters provided in Table 3. Both Henry's law constant ( ) and Langmuir monolayer sorption capacity ( ) tend to increase for highly condensable gases (CO2 > CH4 > N2). We note that the DMS model is highly sensitive to the original pressure range and extrapolated distance, 91 and a detailed analysis of the DMS fitting parameters is beyond the scope of this study. However, as described in      Table 3 according to eqs 5 and 6. As shown in Table 4, the calibrated diffusion coefficients of CO2, CH4, and N2 in PIM-1 are (29.9 ± 4.2) × 10 −8 cm 2 /s, (43.8 ± 9.9) × 10 −8 cm 2 /s, and (184.1 ± 38.1) × 10 −8 cm 2 /s, respectively. Considering the significant scatter in the experimental values, the calibrated diffusion coefficients show reasonable agreement with the experimental values; however, they do not follow the experimental trend 89,90 that CO2 diffusivity is higher than or comparable to N2 diffusivity with CH4 diffusivity being the lowest. Some factors are responsible for these discrepancies. (1) When comparing calibrated diffusion coefficients to experimental values, concentration-independent diffusivity is assumed; however, it does not necessarily hold depending on the gas species. 81 (2) The DMS fitting parameters are highly sensitive to the pressure range, 91 and the values of (eq 7) have significant uncertainties. 85 (3) The TraPPE-UA force field 37,38 used for the nonbonded interaction has been optimized for the quantitative prediction of thermophysical properties, particularly phase equilibria, over a wide range of physical conditions; however this does not guarantee the quantitative prediction of transport coefficients. (4) Although the static structure factor of the simulated PIM-1 are consistent with experimental data ( Figure S1 in Supporting Information), a dynamic chain mobility in the simulated PIM-1 may be different from that in experimental PIM-1 due to a short-chain length of the former (DP = 200) compared to a typical value of the latter (DP > 500); 16,95 this may result in differences of gas diffusion behavior. refers to the diffusion coefficient obtained from MD simulations via eq 4. Each uncertainty represents a standard error obtained from five independent simulations. b refers to the diffusion coefficient calibrated to approximate in eq 5. c The experimental temperature and pressure are 303 K and 0.2 bar, 90 and 308 K and 4 atm. 89 Despite these potential drawbacks, our model systems do provide some qualitative insights into the gas diffusion behavior in the vicinity of PIM-1/silica interfaces. As shown in Table 4 MMM-1 and MMM-2 is due to a strong quadrupole-dipole interaction between CO2 and the hydroxyl groups of the silica surfaces, which retards CO2 diffusion, whereas N2 has a weaker quadrupole moment and CH4 does not. It has been reported that incorporating polar hydroxyl groups into polyimides with intrinsic microporosity improves CO2 adsorption due to their strong dipole moments. 96 The strong quadrupole-dipole interaction between CO2 and surface hydroxyls is also shown in Figure 5, where a significant amount of CO2 is adsorbed even at low pressure in the PIM-1/silica mixing regions. Here, we see two competing effects of silica surfaces: the formation of larger pores, which improves gas diffusion, and the quadrupole-dipole interaction, which retards gas diffusion. More specifically, our findings show that the former has a greater impact on N2 and CH4 diffusivity, whereas the latter has a greater impact on CO2, resulting in a reduction in CO2 diffusivity.
Experimentally, the addition of silica nanoparticles into a PIM-1 matrix leads to an increase in CO2 diffusivity. 20 The difference between our simulation and experimental results is because aggregated silica nanoparticles are dispersed in a polymer matrix with the weight ratio of up to ~ 40 wt% in experiments, 19,20 whereas our MMM models highlight PIM-1/silica interfacial systems. Nonetheless, our findings suggest that the presence of silica surfaces results in different diffusion behaviors of gases in the vicinity of PIM-1/silica interfaces, providing intriguing molecular insights that are unobtainable by experiments. Finally, we opine that to model more realistic MMMs considering the dispersion of aggregated nanoparticles in a PIM-1 matrix, a coarse-grained (CG) MD technique would be an attractive choice because CGMD allows for simulation at a much larger spatiotemporal scale than can be achieved with all-atom MD. [97][98][99] However, the CG force field requires precise parametrization to reproduce static/dynamic properties of interest.

CONCLUSIONS
In this study, we investigate the adsorption and diffusion behaviors of CO2, CH4, and N2 in between gas molecules and silica surface hydroxyl groups, which retards gas diffusion. More specifically, the former significantly affects CH4 and N2 diffusivities because N2 has a weak quadrupole moment and CH4 has no quadrupole moment, whereas CO2 diffusion is more susceptible to the latter due to the strong quadrupole-dipole interaction between CO2 and surface hydroxyls. These findings provide intriguing insights into the behavior of gas diffusion in the immediate vicinity of PIM-1/silica interfaces, which are inaccessible to experiments.
We intend to investigate the effect of silica surface modification on gas adsorption and diffusion behaviors in PIM-1/silica interfacial systems in the future. It has been demonstrated experimentally that surface modification of silica nanoparticles leads to the formation of micro-and mesopores in a polymer phase, resulting in the enhancement of gas diffusivity. 29,100 However, its molecular-level mechanisms still remain elusive, and molecular simulations would provide valuable insights that complement experimental observations.

Supporting Information
Structure factors of PIM-1 membranes; force-field parameters of the Morse-style potential used for constructing amorphous silica structure; structural characteristics of amorphous silica; CWCA force-field parameters for silica surfaces in MMMs; system sizes of model MMMs before/after 21-step equilibration; the relationship between chemical potential and pressure for the bulk gas phases of CO2, CH4, and N2; and gas density distributions in MMM-2 obtained from GCMC simulations (PDF)

Author Contributions
# Y.Y. and Y.T. equally contributed to this work.

Notes
The authors declare no competing financial interest.

ACKNOWLEDGMENT
This work was partly supported by Research and Education Consortium for Innovation of Advanced Integrated Science (CIAiS) and JSPS KAKENHI Grant No. 20K14644.

S1. STRUCTURE FACTORS OF PIM-1 MEMBRANES
The structure factor is related to a radial distribution function as 1 where is the interatomic distance and is the number density. can be measured by neutron or X-ray scattering experiments.   Table S1.    (Tables S2 and S3).

S4. RELATIONSHIP BETWEEN CHEMICAL POTENTIAL AND PRESSURE
The chemical potential can be expressed as follows: where is the ideal gas contribution and is the excess chemical potential due to energetic interactions. The former is given as where is the Boltzmann constant, Λ is the thermal de Broglie wavelength, is the pressure, and is the temperature.
To account for the nonideality of the gases at high pressure, bulk gas-phase GCMC simulations at 300 K in the μVT ensemble are performed to obtain the relationship between the chemical potential and pressure for CO2, CH4, and N2. In each GCMC simulation, the chemical potential ( ) is specified as an input parameter using eq S4 with being the target pressure (1, 2, 5, 10, 20, 30, and 40 atm); in reality, a final pressure obtained from the GCMC simulation is different from the target pressure due to the effect of . To ensure reliable statistics for each gas species, the simulation box size is adjusted for each such that the system contains at least 100 molecules after the GCMC simulation has sufficiently converged. Figure S5 depicts the obtained relationship between the pressure and chemical potential.
CH4 and N2 reasonably follow eq S4 up to a pressure of ~ 40 atm, showing only slight differences between and (i.e., small ). Meanwhile, CO2 begins to deviate from eq S4 with increasing pressure (> 20 atm), showing a nonnegligible effect of . Table S8 summerizes the relationship between the chemical potential and pressure for each gas species. Figure S5. Relationship between the pressure and chemical potential for CO2, CH4, and N2. The solid, dashed-dotted, and dashed lines represent ideal gas relations for CO2, CH4, and N2, respectively, calculated using eq S4.