Relationship between microporous structure and light gas transport through glassy polymeric membranes revealed by molecular simulations

Abstract

Microporous glassy polymers are attractive materials for gas separation membranes, due to their high permeability and tailorable selectivity, provided physical aging can be delayed. The archetypal microporous glassy polymer, PTMSP, can be blended with a hyper-crosslinked isatin–triptycene porous polymer network (PPN) to delay physical aging. However, while PPN is effective at reducing physical aging, it also affects the permeability of light gases through PTMSP. Molecular dynamics simulations were used here to shed fundamental light on the mechanisms responsible for these effects. Atomistic models are developed that satisfactorily reproduce experimental observations such as matrix density and cavity size distributions for neat PTMSP as well as for PTMSP–PPN blends. Analysis of the simulation results suggests that physical aging is delayed because the PPN inclusions slow down PTMSP relaxation while reducing the connectivity between free volume pockets. To understand how PPN inclusions affect light gas permeability, the atomistic models developed are used to probe CO₂ and CH₄ diffusion and sorption. These simulations are conducted for PTMSP matrices exhibiting varying density, towards reproducing experimental permeability data. Interrogating the simulation trajectories, it is discovered that while CH₄ travels preferentially through the free volume cavities, CO₂ preferentially interacts with the available surfaces, especially in the presence of PPN. These differences help interpret experimental observations. The transport of C₂H₆ and H₂S through PTMSP matrices was also investigated. These gases, and in particular H₂S, were found to absorb within PTMSP, yielding very low diffusion coefficients. The ability to predict differences in diffusion pathways for various natural gas components unveils the possibility of engineering membranes to control permeability and selectivity towards large-scale applications.

---

1. Introduction

Industrial separations have been associated with 10–15% of the world energy consumption.¹ It has been estimated that replacing thermal separation processes, i.e., distillation, with membrane-based technologies could reduce the energy intensity by as much as 90%. Achieving this goal would contribute to the United Nation sustainable development goal #7.² Polymeric membranes have attracted attention because they are affordable, can be produced at scale, and can be engineered to achieve desired permeability, selectivity, and stability.³ However, polymeric materials are intrinsically dynamic systems, whose properties can change upon exposure to heat, oxygen, radiation, and UV light. The decay in performance that may arise over time is referred to as aging.^4,5 Physical aging primarily results from an increase in density and a corresponding reduction in free volume;^6–10 chemical aging arises from partial oxidation due to exposure to molecular oxygen;^11–13 blocking aging is due to the uptake of non-volatile environmental impurities, which reduces free volume accessibility.¹³ To appreciate the impact of aging, it helps to remember that the oxygen permeability in poly(trimethylsilylpropyne) (PTMSP) can decrease by up to 2 orders of magnitude when the polymer is stored at ambient conditions for 4 years;⁹ prolonged storage leads to polymer oxidation, which could rupture the molecular chains;¹³ heating PTMSP films in air at 100 °C for 15 hours decreases oxygen permeability by an order of magnitude.¹⁴

PTMSP is the prototypical polymer used in aging studies.¹⁵ This polymer exhibits a glass transition temperature (T_g) above 350 °C and an exceptionally high free volume fraction, up to 30%.^16,17 The latter property, coupled with the polymer rigid backbone, yields extraordinarily high permeability to gases and organic vapors.¹⁸ Its microporosity arises from a unique molecular architecture, which includes a stiff carbon–carbon backbone with alternating single and triple bonds and bulky trimethylsilyl side groups. These structural features introduce substantial steric hindrance, preventing efficient chain packing and resulting in large, interconnected free volume cavities. Strategies that can mitigate aging in high free volume glassy polymers include the incorporation of non-porous and porous nanomaterials,^19–26 post-modification of membranes through polymer cross-linking,^27–29 the incorporation of nanoparticles in the polymer matrix,^24,30 and the integration of various techniques.^30–32

In a prior experimental study,²⁶ our team investigated the properties of polymer matrices composed of PTMSP and a porous polymer network (PPN) based on triptycene and isatin.^33,34 These PPNs were selected due to their high surface area, tailorable microporosity, excellent compatibility with PTMSP, and ability to delay aging. Indeed, over a three-week period, N₂ permeability decreased by 40% in neat PTMSP, while it only decreased by 15% when PTMSP was blended with 5% PPN. Our experimental data showed that the PTMSP microporous structure changed upon PPN addition, which was accompanied by changes in gas permeability and selectivity. Understanding the molecular mechanisms responsible for these observations, and in particular the relationship between microporous polymer structure and gas transport could help design separation membranes with enhanced performance. To achieve this goal, we conduct molecular dynamics simulations for systems of composition similar to those considered experimentally, and we use the experimental data to validate our computational protocols.

Molecular dynamics (MD) has been successful in predicting the sorption and transport rates of gases and solvents in both rubbery and glassy polymers.³⁵ Atomistic models have been constructed for a polymer of intrinsic microporosity, PIM-1,³⁶ and then used to predict the permeation of small gases, N₂ sorption isotherms at 77 K, and free-volume distributions. The transition state theory was used to calculate solubility and diffusion coefficients, obtaining results in good agreement with experiments for most small gases tested, except CO₂. Ozcan et al.³⁷ conducted concentration gradient-driven MD (CGD-MD) simulations to investigate the transport of H₂ and CH₄ through PIM-1/ZIF-8 membranes. Their results, based on an accurate atomistic representation of polymer/metal–organic framework (MOF) composites, serve as a powerful tool for characterizing and understanding gas transport and separation mechanisms in mixed matrix membranes. MD was also employed to study structure and transport in porous materials, where molecular interactions with pore walls and nanoscale confinement control diffusion,^38–41 as well as to investigate the molecular transport mechanism of gases through amorphous polyethylene (PE),⁴² shedding light on the molecular mechanisms responsible for transport. Building on this foundation of successful efforts, once atomistic PTMSP and PTMSP–PPN models are validated against available experimental data, they can be used to predict the behavior of light gases such as CH₄ and CO₂, for which experiments are available, H₂S, a common contaminant in natural gas, as well as other gases, such as C₂H₆. It is particularly attractive to predict sorption and transport of toxic gases, such as H₂S, in membranes without running experiments. The overarching goal of this study is to analyze the molecular trajectories simulated for each of the gases under different conditions, to discover viable mechanisms that control permeability and selectivity in separation systems.

The remainder of the manuscript is organized as follows: in Section 2, we describe the simulation methods and algorithms; in Section 3, we present a summary of the results we obtained, and we discuss their relevance in relation to experimental observables. In Section 4, we conclude with an overview of our observations, their implications, and potential suggestions for future research.

2. Methodology

The simulation methodology is discussed in detail in the ESI.† In short, the initial steps involved (a) the construction of one short PTMSP polymer chain comprising 101 repeat units, with cis- and trans-arrangement of the monomers distributed randomly along the chain (cf. Fig. 1a), and (b) the creation of one porous polymer network (PPN) obtained by incorporating 22 units of triptycene and 45 units of isatin distributed randomly along the chain (cf. Fig. 1b). The components just described were used to build two distinct systems to explore the structural characteristics of the polymeric membranes. The first system, ‘neat PTMSP’, consists of 20 PTMSP polymer chains; the second system, ‘PTMSP–20PPN’, contains 20 PTMSP chains and 5 PPN molecules. Higher PPN content was considered not relevant to compare against the experimental system considered in our prior work,²⁶ while lower PPN content was considered not to be sufficient to provide statistical accuracy. To ensure that the results are robust, five independent initial configurations were prepared for each system. This protocol was necessary to reasonably capture the variability in molecular structure expected across polymeric membranes.

Fig. 1 Atomistic models for (a) one PTMSP polymer chain consisting of 101 monomer units and (b) one PPN nanoparticle composed of 22 triptycene units and 45 isatin units. The color scheme is as follows: carbon (cyan), hydrogen (white), oxygen (red), nitrogen (yellow), and silicon (orange). The molecular structures of the repeat units are shown for completeness.

The simulation protocol was meant to replicate the solution-casting technique implemented experimentally.²⁶ Both PTMSP and PPN were dissolved in toluene, which was then gradually removed as the simulations were conducted in the NPT ensemble (constant number of particles, pressure, and temperature). The algorithm, described in Section S1 of the ESI† was implemented to produce both PTMSP and PTMSP–20PPN membranes. Fig. 2 presents snapshots of the final configurations for both matrices, while Table 1 provides the cubic simulation box size and the density of the simulated systems. When interpreting the simulation results, one must always consider possible computational limitations. For example, the size of a simulation box is limited by the computational resources available. To replicate the properties of polymeric membranes, it is necessary that the simulation box is larger than the size of connected free volume regions within the material, to prevent the availability of unphysical preferential transport pathways.

Fig. 2 Simulation snapshots illustrating the final configurations of (a) PTMSP and (b) PTMSP–20PPN matrices. The color scheme is as follows: cyan (carbon), white (hydrogen), and orange (silicon). For improved visual distinction, atoms in the PPN are colored purple. The snapshots reveal amorphous morphology and the presence of voids within both model membrane structures.

Table 1 Simulated densities of the PTMSP and PTMSP–20PPN polymer matrices (the densities were calculated by averaging the results over five independent samples for each system). Experimental data are from ref. 26. The simulations are conducted in cubic boxes

Material	Box length (Å)	Density (g cm^−3)
		MD	Experiment
PTMSP	79.28 ± 0.08	0.758 ± 0.002	0.759 ± 0.004
PTMSP–20PPN	81.82 ± 0.40	0.866 ± 0.010	0.864 ± 0.012

---

Once the simulated density was comparable to experimental values for fresh membrane samples, subsequent simulations were performed in the NVT ensemble (constant number of particles, volume, and temperature – i.e., the density remained constant) at 298 K, enabling the characterization of the porous structure. Additional simulations were performed in the NPT ensemble (constant number of particles, pressure, and temperature) at various temperatures to comprehensively examine the aging behavior, starting from the matrices described in Table 1.

Forcefields and algorithms implemented to produce the simulated trajectories as well as those used to analyze the simulation results are discussed in ESI.†

3. Results and discussion

3.1 Reproduction of experimental properties

The results presented in this section assess whether the models developed for PTMSP and PTMSP–PPN systems reproduce the experimental data reported previously for these materials.²⁶ When the agreement is satisfactory, the simulated trajectories are interrogated further to extract additional insights, which complement experimental observations, and provide molecular-level interpretation for the experimental results.

3.1.1 Polymer matrix structure. The cavity size distributions (CSDs) obtained from the simulations of PTMSP and PTMSP–20PPN matrices are shown in Fig. 3. To complete the comparison against experimental properties, we also computed the CSDs for a system composed only of PPN. As discussed in ESI,† the average density of PPN systems used for these calculations was 1.21 ± 0.04 (g cm⁻³). The CSDs obtained exhibit multiple peaks, indicating a range of cavity sizes, with relatively broad distributions, spanning from 0.5 to 3 nm, with prominent peaks at 0.5, and between 2 and 3 nm. The CSD results for PPN show a significant peak for sizes less than 0.5 nm. Both PPN and the two membranes were experimentally tested using Positron Annihilation Lifetime Spectroscopy (PALS); the results, reported in our previous study,²⁶ showed distinct peaks at ∼0.55, 1.3, and 1.74 nm, estimated by fitting the experimental annihilation lifetimes of positrons with the empirical Tao–Eldrup model.^43,44 Considering the differences between the two approaches, the data variability due to the small size of the simulated system (cf. the variability of the profiles obtained for the five different samples considered, shown in Fig. S1 of ESI†), the uncertainties associated with both experiments and simulations, and the fact that the densities of simulated vs. experimental samples are similar (cf. Table 1), it is concluded that the simulation models provide a reasonable representation of the experimental systems because the cavity sizes observed in the simulations are in the same order of magnitude as those extracted by fitting PALS data. An important observation is that the incorporation of 20 wt% PPN into PTMSP shifts the cavity size distribution towards smaller cavities, reducing the prevalence of cavities in the 2–3 nm range and somewhat enhancing the formation of smaller-sized cavities (less than 1 nm). It is noted that the PALS data suggest small cavities (less than 0.5 nm and of ∼1 nm in size) for the PPN inclusions, but do not support changes in the prevalence of cavities smaller than ∼1 nm upon blending PPN with PTMSP. This difference could be important in controlling the transport of lights gases, especially those of small size, although it does not seem to have affected the results obtained here, as discussed later.

Fig. 3 Cavity size distributions observed from simulations of (a) PPN, (b) PTMSP and (c) PTMSP–20PPN matrices. The distribution in panel (a) represents the averaged data from three independent runs. The results in panel (b) and (c) represent averages over five independent samples. Detailed results for the five replicates of each simulated system are provided in Fig. S1 of the ESI.† Corresponding experimental data are reported in ref. 26.

Both the interpretation of PALS experiments and the CSDs in Fig. 3 consider spherical cavities. However, snapshots (e.g., cf. Fig. 2) suggest that the cavities in the simulated systems are not spherical. Hence, we complemented the CSDs results with constrained stochastic ray tracing analysis. Ray tracing analysis is obtained by quantifying the percent of ‘rays’ of a given length that can be inserted in the simulated system without overlapping with the polymer/PPN molecules. If the cavities in the system were spherical, results from ray tracing analysis would be similar to those obtained using CSDs. If the cavities were not spherical, the CSDs would identify the spheres with higher diameter that can be placed in the system without overlapping with the polymer, while ray tracing would capture the possibility of inserting rays that span the longest diameter within a cavity. Hence, this method provides a statistical representation of the free volume topology, capturing details of cavity shape, surface texture, and spatial connectivity. These results help interpret the degree of cavity accessibility and continuity, which are key factors influencing gas transport through polymeric matrices. The results are presented in Fig. 4 for both PTMSP and PTMSP–20PPN. While the results in Fig. 3 do not show cavities of size larger than ∼3.3 nm, the distributions in Fig. 4 show that rays of length 4 nm or longer (up to 8 nm for the PTMSP sample) can be observed in the simulated samples. The introduction of PPN yields a shift toward shorter ray lengths, accompanied by the decrease of the shoulder at ∼4 nm, and the enhancement of the peak at ∼0.6 nm. The combination of Fig. 3 and 4 suggests that adding 20% PPN yields a reduction in cavity connectivity compared to neat PTMSP. The constrained ray tracing distribution for PTMSP–20PPN exhibits a characteristic asymmetric peak, contrasting with the broader profile of neat PTMSP. This distinct peak is consistent with a structural rearrangement in the cavity network, favoring smaller, less interconnected voids upon the addition of PPN to PTMSP.

Fig. 4 Constrained stochastic ray trace distribution results as a function of ray length obtained for PTMSP and PTMSP–20PPN. The results were averaged over five independent samples for each simulated system.

To appreciate the structural differences, we present in Fig. 5 three-dimensional visualizations of accessible cavity space distribution and arrangement within both polymer matrices. In PTMSP (cf. Fig. 5a), a highly interconnected network of micropores, characterized by rather expansive accessible regions, is visible. These structures are consistent with the low density of the material and could provide enhanced diffusion pathways for gases, consistent with the high gas diffusion coefficients reported experimentally for freshly cast PTMSP samples, relative to PTMSP–20PPN. Conversely, PTMSP–20PPN (cf. Fig. 5b) exhibits a microporous structure less interconnected compared to that found in PTMPS. In Fig. 5c, we show the distribution of PPN inclusions within the matrix. Comparing qualitatively the cavities distribution in PTMSP–20PPN with the accumulation of PPN inclusions, it appears that the PPN effectively ‘choke’ the connectivity between the PTSMP micropores. This qualitative visualization supports the quantitative findings from both CSD (cf. Fig. 3) and stochastic ray tracing analyses (cf. Fig. 4).

Fig. 5 Visualization of micropore connectivity for (a) PTMSP and (b) PTMSP–20PPN matrices. The surface of the void spaces is represented in blue, highlighting accessible cavities regions within the polymer matrix, while the volume occupied by atoms of polymer and PPN nanoparticles are colored in yellow. The surface mesh was calculated by the alpha-shape method⁴⁵ with a probe sphere radius of 0.3 nm as available in the software Ovito.⁴⁶ To highlight the effect of PPN nanoparticles, panel (c) illustrates the molecular structure of PTMSP–20PPN, with purple atoms representing the PPN. The color scheme for PTMSP is that used in Fig. 2.

3.1.2 Dynamical properties of the polymeric matrices. Although aging at ambient conditions occurs experimentally over a period of at least hours or days, atomistic MD can reveal useful information regarding the effect of PPN inclusions in delaying the mobility of PTMSP chains. It should be recognized that the equations of motion in atomistic MD simulations cannot at present be integrated for times equivalent to hours and days for large systems. However, the aging process can be accelerated by conducting simulations at temperatures higher than experiments (this would correspond to thermal aging), and by considering shorter polymers than those used experimentally, thus reducing the effective viscosity of the systems. We interrogate the simulated trajectories to (a) analyze atomic-scale dynamics of the polymer chains via velocity–velocity autocorrelation functions (VACFs), (b) quantify the fractional free volume (FFV) as a function of temperature, and (c) analyze the atomic details of PTMSP–PPN interactions.

In Fig. 6, we present VACFs obtained for three distinct atomic groups, chosen to capture a range of dynamic behaviors within the PTMSP matrix: (A) the carbon atoms in the CH₃ groups of trimethylsilyl groups, (B) the carbon atoms in the pendant CH₃ groups, and (C) the carbon atoms in the polymer backbone. The C–Si bonds show the most restricted dynamics, as indicated by the rapid VACF decay in Fig. 6a, reflecting limited mobility. The pendant groups (Fig. 6b) exhibit intermediate dynamics, characterized by an oscillatory behavior that suggests sensitivity to local interactions. The backbone carbons (Fig. 6c) exhibit persistent oscillatory behavior. Notably, the VACF profiles do not change upon incorporation of PPN, across all groups analyzed. These results suggest that while PPN alters the matrix morphology and reduces cavity connectivity, it does not impede the intrinsic segmental mobility of the polymer, thereby maintaining local vibrational dynamics critical for membrane performance. The simulated VACF results are consistent with prior FT-IR experimental data,²⁶ which did not show changes in the spectra for C=C and Si–CH₃ bonds upon adding PPN to the blends, qualitatively consistent with our simulation results. This observation, which is consistent also with longer-time results experimentally probed by the NMR data reported in our prior study,²⁶ underscores the importance of further investigating the role of PPN in mitigating the aging process in PTMSP membranes.

Fig. 6 Velocity–velocity autocorrelation functions (VACFs) for carbon atoms in (a) C–Si bonds, (b) pendant groups, and (c) polymer backbone in PTMSP and PTMSP–20PPN. The results for the different materials (red vs. black) are indistinguishable from each other. The results were averaged over five independent samples for each simulated system. These simulated results are consistent with experimental FTIR data reported in ref. 26.

In Fig. 7, we present simulated fractional free volume (FFV) profiles as a function of probe radius and temperature. The simulations are conducted in the NPT ensemble, where the system pressure is maintained constant. As the polymers density increases, the simulation box size decreases and the FFV decreases as well, mimicking the process of physical aging. The results shown in Fig. 7 are obtained by initiating the simulations with the configurations used to achieve the densities shown in Table 1. The temperature is set at the desired conditions, and the FFV is monitored as the simulations progress. As mentioned above, conducting these simulations at elevated temperatures helps speed up aging, which resembles thermal densification and structural relaxation. Each curve is the average of five independent simulations. FFV profiles reflect the proportion of unoccupied or accessible space within the polymer matrix relative to the total volume of the sample.⁴⁷ As the probe radii (R_p) increases from 0 to 0.25 nm, the FFV for both membranes decreases, indicating that larger probes face greater difficulty in accessing the small voids present throughout the polymeric matrix. This trend is consistent across all temperatures, although higher temperatures yield lower overall FFV profiles due to polymer matrix compaction. The observed decrease in FFV with increasing temperature is associated with physical aging. Indeed, at elevated temperatures, physical aging is accelerated, causing the equilibrium state of the material to shrink.⁴⁸ It is noted that the simulated changes in FFV occur much faster than experimental physical aging of PTMSP samples. It is hypothesized that this discrepancy is a consequence of the low molecular weight of the PTMSP chains simulated, which are too short to yield physical phenomena such as entanglements which will certainly delay aging. This hypothesis, however, cannot be tested with current computational resources.

Fig. 7 Fractional free volume for PTMSP and PTMSP–20PPN matrices as a function of probe radius (panels a & b) and temperature (panel c). The results in panel c are obtained for probes of size 0.043 nm. These simulation results are qualitatively consistent with prior experimental data.²⁶

Consistent with literature,⁴⁹ the results in Fig. 7 depend on the probe size used for sampling the simulated systems. However, the results are consistent across all probe sizes, and aligns with our experimental findings.²⁶ It has been recommended by Hofmann et al.⁵⁰ to use results from probes of size 0.043 nm to compare FFV from simulation to PALS experiments. We use this probe size to quantify the effects of PPN on aging via detailed analysis of Fig. 7 results. At the limiting probe size of 0.043 nm, the FFV of PTMSP at 298 K starts from ∼0.357 ± 0.002 while that for PTMSP–20PPN is ∼0.289 ± 0.008. It is recognized that the experimentally measured decrease upon incorporating PPN (∼15.8%) is slightly lower than that observed in simulations for probes of size 0.043 nm (∼18.9%). As discussed in what follows, a discrepancy between experiments and simulations is expected. Overall, the decrease in FFV suggests that the addition of PPN leads to a tighter packing within the polymer matrix. As the temperature increases from 298 K to 598 K, the FFV of the neat PTMSP decreases from 0.357 ± 0.002 to 0.269 ± 0.016, indicating significant thermal densification. In contrast, the reduction in FFV for PTMSP–20PPN is less pronounced, demonstrating that the PPN mitigates the loss of free volume upon aging. Although the absolute FFV obtained for PTMSP–20PPN is lower than that obtained for PTMSP, the smaller temperature-dependent variation in FFV is interpreted as enhanced structural stability.

In Fig. 8, we provide simulation snapshots suggesting that PPN stabilizes the PTMSP chains by interpenetrating and interacting with the polymer via small surface protrusions on the PPN inclusions. Overall, the results presented so far suggest that PPNs do not affect the short-term relaxation properties of the PTMSP chains (cf., Fig. 6). However, the PPN inclusions have a pronounced effect on reducing the FFV within the matrices, an effect that is more pronounced at low temperatures (cf., Fig. 7). It is useful to point out that as temperature increases, the FFV in PTMSP decreases monotonically, while that in PTMSP–20PPN decreases only slightly. This suggests that PPN inclusions slow down PTMSP chains dynamics, which is consistent with enhanced aging resistance. The snapshot in Fig. 8 suggests that this reduced polymer mobility might be due to increased friction between inclusions and polymer chains, as PPN inclusions seem to clamp the PTMSP chains in between them. If proven, this hypothesis could be used to engineer new composite PTMSP-based membranes with enhanced resistance to physical aging.

Fig. 8 Snapshots representing PTMSP–PPN interactions: (a) one PTMSP chain (violet) is captured between two PPN inclusions (green and red); (b) details of the interface between one PPN inclusion (orange) and one PTMSP chain (cyan).

3.2 Prediction of light gases behavior

The results presented in this section leverage the observation that the simulation models developed here reproduce reasonably well the properties of the experimental polymeric matrices; hence, atomistic molecular dynamics simulations were used to predict thermodynamic and transport behavior of some light gases. The gases considered, C₂H₆, CH₄, CO₂, and H₂S, are common in natural gas. Prior experimental data are available for CH₄ and CO₂ in PTMSP and PTMSP–PPN matrices, thereby allowing us to further validate the approach. For the other two gases, H₂S being a toxic contaminant, the simulations are used as predictive tools. The reliability of the predictions relies on the good agreement with experimental data obtained for CH₄ and CO₂.

3.2.1 Solubility of light gases in the polymeric membranes. Using the PTMSP and PTMSP–20PPN matrices representing freshly cast materials (i.e., the simulation boxes used to start the investigation of FFV changes with temperature shown in Fig. 7), it is possible to compute the partition coefficients of C₂H₆, CH₄, CO₂ and H₂S between polymer membrane and gas phase at infinitely low pressure at 298 K (i.e., infinite dilution sorption coefficients, cf. Fig. 9a). These simulations are conducted at constant volume, so that the polymer matrices do not compact during the simulations. Since these calculations are conducted at low gas density, this approximation is considered reasonable. Should high gas densities be considered, it is possible that the polymer matrix relaxes and changes its structure due to interactions with the gases.

Fig. 9 Panel (a) represents the simulated partition coefficients of CH₄, CO₂, H₂S, and C₂H₆ between polymeric matrices and the gas phase (at infinitely low pressure). Panel (b) represents the diffusion coefficients simulated for CH₄, CO₂, H₂S, and C₂H₆ in pristine matrices. Results are obtained in neat PTMSP and PTMSP–20PPN systems. All simulations are conducted at 298 K. Panel (c) represents diffusion coefficients predicted for CH₄, CO₂, H₂S, and C₂H₆ in densified PTMSP samples. The simulations conducted for obtaining results shown in panel c were conducted by randomly placing 40 molecules of each gas into densified PTMSP matrices and applying periodic boundary conditions in all three dimensions. The systems were equilibrated in the NVT ensemble at 298 K and 1 bar for 40 ns, followed by a 4 ns production run for diffusivity analysis. The densified PTSMP density is 0.922 ± 0.002 (g cm⁻³).

Simulations show that the partition coefficient of CH₄ in neat PTMSP, ∼15, is the lowest out of all the gases considered. The highest partition coefficient is the one obtained for C₂H₆, which is consistent with ethane having higher critical temperature than methane, as gas sorption in polymers increases with gas condensability.⁵¹ Upon incorporating 20 wt% PPN, the partition coefficient for CH₄ increases to ∼18; similarly, the partition coefficient for C₂H₆ increases from ∼275 to ∼360 upon adding PPN to PTMSP, indicating that PPN inclusions have a relatively small effect on CH₄ and C₂H₆ solubility; however, the partition coefficient for CO₂ increases from ∼35 for neat PTMSP to 78 for PTMSP–20PPN, showing that inclusion of PPN significantly enhances CO₂ sorption. The results obtained for CH₄ and CO₂ exhibit reasonable agreement with experimental data,²⁶ where infinite dilution sorption coefficients increased by 114% for CO₂ and ∼50% for CH₄ when comparing PTMSP–20PPN to neat PTMSP data. In our prior work,²⁶ it was hypothesized that the substantial increase in CO₂ solubility might be due to favorable interactions between CO₂ and the lactam group on PPN. The simulation results show that the increased solubility upon adding PPN to PTMSP is even more pronounced for H₂S, in which case the partition coefficient increases from ∼138 to ∼540. While H₂S is more condensable than CO₂, and, therefore, more soluble in polymers, the simulation results indicate that H₂S is strongly attracted to both PTMSP and PPN chains. Specifically, the larger H₂S condensability (T_c = 100 °C) relative to CO₂ (T_c = 31 °C) justifies the larger partition coefficient in PTMSP; in addition, H₂S is expected to favorably interact with lactam groups in PPN, similarly to CO₂.

To test our hypotheses, we analyze simulation snapshots. In Fig. 10, the four snapshots illustrate the distribution of CH₂, C₂H₆, CO₂, and H₂S molecules throughout the PTMSP–20PPN matrices. The results show that CH₄ molecules can be most likely found within the free volume pockets within the polymer matrices, while C₂H₆ molecules are both in the free volume pockets and absorbed in the polymer. CO₂ molecules have a preference to reside near the surface of the PTMSP–20PPN materials, with preferential interactions with the PPN nitrogen atoms. Different than all other gases, H₂S is found to absorb throughout the PTMSP polymer mass, as well as within the PPN inclusions. To further support our interpretation, we calculate density profiles for C₂H₆, CH₄, CO₂, H₂S, and the PPN nitrogen atoms along planes that cut the PTMSP–20PPN matrices. The results, shown in Fig. S3 of ESI,† reveal a high local concentration of CO₂ and H₂S near the PPN nitrogen atoms (i.e., near the lactam groups), confirming their role in enhancing the solubility of these gases. In contrast, C₂H₆ and CH₄ are more uniformly dispersed throughout the PTMSP–20PPN matrices. This uniform distribution suggests that C₂H₆ and CH₄ solubility is primarily influenced by the physical free volume across the polymer rather than by specific chemical interactions.

Fig. 10 Representative simulation snapshots showing the distribution of (a) C₂H₆, (b) CH₄, (c) CO₂, and (d) H₂S molecules within PTMSP–20PPN matrices. The color scheme for PTMSP and PPN corresponds to that used in Fig. 2. C₂H₆, CH₄, CO₂, and H₂S are represented by blue, yellow, black, and green spheres, respectively.

3.2.2 Transport of light gases in the polymeric membranes. The diffusion coefficients obtained for pure light gases in pristine PTMSP and PTMSP–20PPN are presented in Fig. 9b. The structure of the polymer matrices used for these calculations is shown in Fig. 2. The results show that in PTMSP, CH₄ has the highest diffusion coefficient, slightly above that of CO₂, while H₂S has the lowest, less than half that of CO₂; in contrast, C₂H₆ has diffusion coefficient slightly higher than H₂S. These data are at first sight surprising, because the diffusion coefficients of small molecules in polymers generally decrease with increasing penetrant size. Therefore, CO₂ (kinetic diameter = 3.3 Å) is expected to diffuse faster than CH₄ (kinetic diameter = 3.8 Å). However, in high free volume glassy polymers, such as PTMSP, CO₂ and CH₄ exhibit almost similar diffusion coefficients, as confirmed by both computational and experimental data. Hofmann et al.,⁵⁰ e.g., reported simulated diffusion coefficients of 0.8 × 10⁻⁹ m² s⁻¹ for CO₂ and 1.7 × 10⁻⁹ m² s⁻¹ for CH₄ in pure PTMSP. Similarly, Zhou et al.⁵² reported diffusion coefficients of 3.3 × 10⁻⁹ m² s⁻¹ for CO₂ and 3.8 × 10⁻⁹ m² s⁻¹ for CH₄. Heuchel et al.³⁶ reported calculated diffusion coefficients for CH₄ that exceed those predicted for CO₂ by a factor of ∼35. Experimental studies from Merkel et al. reported concentration-averaged diffusion coefficients of 3.3 × 10⁻⁹ and 3.6 × 10⁻⁹ m² s⁻¹ for CO₂ and CH₄, respectively in neat PTMSP,⁵³ in good agreement with our recently published experimental transport diffusion coefficients at 35 °C and 1 atm (i.e., (3.7 ± 0.3) × 10⁻⁹ m² s⁻¹ for CO₂ and (3.4 ± 0.3) × 10⁻⁹ m² s⁻¹ for CH₄). Considering that the diffusion coefficients in complex matrices depend on the presence of, and access to cavities, the connectivity among cavities, and on preferential interactions between the matrix and the gases,⁵⁴ among other factors, and noting that the polymer microstructure depends on composition and history (cf. Fig. 7), and that the accuracy of describing preferential gas–matrix interactions depends on the force fields, the results obtained here appear to be reasonable. To further test the reliability of the simulated diffusion coefficients, we also simulated an aged (i.e., densified) PTMSP material. This PTMSP sample shows a density of ∼0.922 ± 0.002 (g cm⁻³) (cf. Table S1†), within which only free volume pockets of size in the range 0.5–1.0 nm can be observed (cf. Fig. S2†). In this densified PTMSP model, our simulations (Fig. 9, panel c) predict that CO₂ has the highest diffusion coefficient, followed by CH₄, H₂S, and C₂H₆. These data agree with our prior experiments for CH₄ and CO₂, as well as with expectations based on differences in kinetic diameters between the four gases.

The results shown in Fig. 9, panel b, predict a considerable reduction in the diffusion coefficients of all gases upon the addition of 20 wt% PPN. The results are quantified in a reduction of ∼76% for CH₄, ∼63% for CO₂, ∼61% for H₂S, and ∼55% for C₂H₆, suggesting that the incorporation of PPN leads to an increase in tortuosity and reduce the free volume within the membrane. This result is consistent with the structural results presented above (cf. Fig. 5). It is noteworthy that the reduction in CO₂ diffusivity observed in our simulations is in agreement with experiments, in which transport CO₂ diffusivity decreased by 52% upon incorporation of 20 wt% of PPN in PTMSP.²⁶ The combination of all diffusion coefficients predictions in Fig. 9 highlights that, although PPN incorporation leads to a moderate reduction in the initial free volume, PTMSP–20PPN promises to maintain gas transport properties over time compared to densified PTMSP, supporting the strategy of using PPN inclusions to enhance long-term membrane performance. Moreover, the decrease in diffusivity due to PPN incorporation is compensated by the sorption coefficient increase, which, based on the solution-diffusion model, enhances permeability.²⁶

To gain further insights, we analyzed representative trajectories of the pure gas molecules (cf. Fig. 11). The results reveal distinct diffusion patterns. CH₄ molecules preferentially occupy micropores, within which their trajectories exhibit uninterrupted movement, particularly in PTMSP. This behavior is diametrically opposite compared to that of H₂S, which is absorbed by the PTMSP matrix and therefore diffuses slowly even in pristine microporous PTMSP. Our results show that C₂H₆ has a behavior somewhat in between that of CH₄ and that of H₂S, exhibiting strong attraction with the PTMSP polymer matrices while also diffusing in the free volume cavities, while CO₂ exhibits more constrained diffusion, characterized by shorter and more localized pathways, reflecting stronger interactions with the polymer surface and the PPN inclusions. The trajectories show that PPN affects the diffusion mechanism of CH₄ and C₂H₆ by reducing the connectivity between micropores, and that of CO₂ and H₂S by introducing preferential adsorption sites that delay diffusion. To quantify the effect of temperature on the predicted diffusion coefficients, we conducted additional simulations at 318 K and 348 K, with results presented in Fig. S4 of the ESI.† The diffusion coefficient in neat PTMSP exhibits a pronounced increase with temperature for all gases. Gas diffusion in PTMSP–20PPN is systematically slower than in neat PTMSP, and exhibits a much weaker temperature dependence. Within the small temperature range tested, the temperature dependence of CO₂ diffusion is evident in neat PTMSP, but in PTMSP–20PPN, the increase is much less pronounced, suggesting that these small changes in temperature are not sufficient to affect the preferential PPN–CO₂ interactions. The temperature dependence of H₂S diffusion is less pronounced than that of the other gases in PTMSP, suggesting that its mobility is more constrained by the rigid, microporous structure of the polymer matrix.

Fig. 11 Molecular trajectories of CH₄ (a) and (f), CO₂ (b) and (g), H₂S (c) and (h), and C₂H₆ (d) and (i) within the neat PTMSP membrane (e) and the PTMSP–20PPN blend (j) over a 4 ns production run. The trajectories are color-coded to represent the time evolution of the gas molecules. Panels (e) and (j) highlight the void spaces within the membranes where the trajectories were monitored, with void regions represented in gray, the PTMSP polymer matrix in yellow, and PPN in red for clarity. All visualizations were generated using the software Ovito.⁴⁶ Simulations conducted at 298 K.

Our simulations predict fast CH₄ diffusion in PTMSP (cf. Fig. 9, panel a) because the free volume cavities span through the simulation boxes, leading effectively to infinitely long interconnected cavities. These cavities are highly unlikely to be found in real polymeric materials, and a much larger simulation box than those considered here would be necessary to model realistic scenarios. Since polymers such as PTMSP and PTMSP–20PPN are heterogeneous, as shown in Fig. 2, with cavities separated by dense polymer regions, one could extract diffusion coefficients in the cavities as well as in the dense polymer matrix, and then use analytical or stochastic models to predict the effective diffusivity.⁵⁵ Those approaches are not attempted here because a full three-dimensional description of the matrices is not available.

3.2.3 Preferential gas–polymer matrix interactions. To further probe the molecular mechanisms that control the behavior of light gases in PTMSP membranes, residence autocorrelation functions, C_R(t) were calculated.^38,56,57 The results provide an estimate of how long a molecule remains at a specific location from a reference site. In general, the faster C_R(t) decays from 1 to 0, the quicker the probe molecules leave regions of interest. Our analysis focuses on C₂H₆, CH₄, CO₂, and H₂S molecules residing within a cutoff distance of ∼1.0 nm from the nitrogen atoms of PPN. The cutoff distances were determined based on radial distribution function (RDF) calculations. The results for both RDF and C_R are presented in Fig. 12. C₂H₆ and CH₄ molecules are not strongly attracted to the N atoms in the PPN lactam groups, as shown by the RDF results, and exhibit rapidly decaying C_R(t), indicating that these molecules spend short times near the N atoms of PPN, as expected from Fig. 10. In contrast, CO₂ molecules are strongly attracted to the N atoms, as shown by the pronounced peaks in the correspondent RDF data, and their C_R manifest a slower decay, suggesting prolonged residence times, which contributes to its reduced diffusion. Interestingly, although H₂S shows the highest first peak in the RDF profiles – indicating a strong initial attraction to nitrogen atoms – its C_R(t) decays more quickly than that of CO₂. This suggests that while H₂S interacts strongly with PPN, these interactions are more transient, resulting in shorter residence times compared to CO₂. It is possible that strong H₂S–H₂S preferential interactions effectively reduce the lifetime of H₂S–PPN association, but this possibility has not been tested in detail.

Fig. 12 (a) RDF between gas molecules and N atoms in PPN and (b) residence autocorrelation functions C_R(t) for gas molecules within cutoff distances of ∼1.0 nm from N atoms in the polymer blends. The simulation results were obtained for pure gases in PTMSP–20PPN blends at 298 K, with data averaged over five independent samples.

Finally, vibrational spectra of CO₂ molecules in the polymer matrices were predicted. The results were differentiated based on proximity to N atoms of PPN, as shown in Fig. 13. The power spectra were obtained via a Fourier transform of velocity autocorrelation functions.⁵⁸ For these calculations, atomic velocities were recorded every 4 fs over a total simulation time of 40 ps. The spectral analysis reveals that while the fundamental vibrational modes of CO₂ remain unchanged compared to those observed in the bulk phase,⁵⁹ the intensity differences between spectra suggest confinement-induced variations in molecular behavior. In particular, the asymmetric mode becomes significantly more pronounced when CO₂ is within 1 nm from N atoms, because of preferential interactions. These findings indicate that the incorporation of PPN alters CO₂ dynamics, impacting its diffusion and solubility within the polymer matrix. Experimentally, it is expected that the overall FTIR signal would change with CO₂ loading, as the first molecules would preferentially interact with the N sites; once these are occupied, additional CO₂ molecules would be found throughout the matrix.

Fig. 13 Power spectra for CO₂ molecules within the polymer matrix: (a) within a 1 nm cutoff distance from N atoms and (b) beyond 1.3 nm from N atoms. These cutoff distances were determined based on RDF calculations, as shown in Fig. 12. The spectra were derived from molecular dynamics simulations of CO₂ in PTMSP–20PPN blends at 298 K and 1 bar and averaged over five samples. Shaded regions represent standard deviation.

4. Conclusions and outlook

This study provides molecular insights into the relation between the microstructure of microporous glassy polymers and the transport behavior of light gases. Specifically, both neat PTMSP and PTMSP blended with a triptycene–isatin porous polymer network (PPN) were considered as model materials. Using molecular dynamics simulations, an atomistic model was derived for both polymeric membranes. The results show that incorporating 20 wt% PPN into PTMSP significantly alters the cavities morphology and reduces the fractional free volume. Cavity size distribution (CSD) analysis demonstrates that the addition of PPN shifts the cavity distribution towards smaller sizes, reducing the prevalence of larger cavities (2–3 nm) and increasing the concentration of smaller cavities (<1 nm). This structural transformation occurs while also disrupting the connectivity between larger cavities and enhancing the packing density of the polymer matrix. These results are consistent with experimental data,²⁶ which suggested PPN effectively reduces physical aging of PTMSP. These structural properties reproduce experimental results reported previously,⁶⁰ including those recently published by our group.²⁶

Building on the reliability of the atomistic model, molecular dynamics simulations were used to predict the behavior of light gases, representing common components, including contaminants, of natural gas: ethane, methane, CO₂ and H₂S. Each gas shows unique behavior, dictated by its molecular size, specific interactions with PTSMP and/or PPN, and availability of free volume cavities. Results obtained for CH₄ and CO₂ show marked differences, which are in general agreement with prior experimental data. Predictions for C₂H₆ and H₂S demonstrate the ability of molecular simulations to interrogate systems that can be difficult or risky to handle in the laboratory. The differential diffusion pathways observed for different gases in polymer membranes are fully compatible with the solution-diffusion mechanism.

Our findings highlight an important trade-off in membrane design: while the incorporation of PPN reduces the initial fractional free volume and gas diffusivity due to the disruption of large free volume cavities, it effectively stabilizes the polymer matrix against long-term densification. Moreover, PPN introduces preferential adsorption sites, which may enhance gas selectivity by promoting stronger interactions with specific gas molecules. As a result, the PTMSP–20PPN blends can achieve a balance between permeability and durability.

While the simulation predictions always depend on the accuracy of forcefield models, faithful description of the polymeric matrix, and adequate sampling of the properties of interest, the results presented here suggest that if the three-dimensional matrix structure was known with high accuracy, analytical and/or stochastic methods such as the effective medium theory or kinetic Monte Carlo approaches could be used to predict gas permeability in polymeric mixed matrix materials leveraging molecular dynamics simulation results as input parameters. The conceptual workflow presented, and the results obtained help discover new design strategies to engineer polymeric materials with desired morphology and cavity size distribution and connectivity, and exhibiting both high performance and durability.

Data availability

Input files, including force fields and simulated trajectories, will be deposited on an institutional repository. This information will be available to the interested readers upon request.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

Financial support was generously provided by the US Department of Energy (DOE, Office of Basic Energy Sciences) under grant no. DE-SC0024554. AS acknowledges financial support from the Asahi Glass Chair in Chemical Engineering at the University of Oklahoma. Generous allocations of computing time were provided by the University of Oklahoma Supercomputing Center for Education and Research (OSCER). We sincerely acknowledge the University of Oklahoma's Office of the Vice President for Research and Partnerships (OVPRP) and the Data Institute for Societal Challenges (DISC) for providing GPU resources. We also extend our gratitude to the Stampede3 at Texas Advanced Computing Center, through allocation TG-CHM240008 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. Additionally, we acknowledge the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/T022213/1, EP/W032260/1 and EP/P020194/1), as well as the High-Performance Computing Center at Oklahoma State University. Finally, the authors are thankful to Dr Thang Ha Nguyen, Research Computing Facilitator at the University of Oklahoma for his technical support.

References

1. D. S. Sholl and R. P. Lively, Nature, 2016, 532, 435–437.
2. IEA, IRENA, UNSD, World Bank, WHO, Tracking SDG 7: The Energy Progress Report, World Bank, Washington, DC, 2024, © World Bank, License: Creative Commons Attribution—NonCommercial 3.0 IGO (CC BY-NC 3.0 IGO).
3. L. M. Robeson, J. Membr. Sci., 2008, 320, 390–400.
4. M. Robert, A. El Kaddouri, J.-C. Perrin, K. Mozet, M. Daoudi, J. Dillet, J.-Y. Morel, S. André and O. Lottin, J. Power Sources, 2020, 476, 228662.
5. F. Zhang, R. Yang and D. Lu, Polymers, 2023, 15, 1928.
6. S. Tasaka, N. Inagaki and M. Igawa, J. Polym. Sci., Part B: Polym. Phys., 1991, 29, 691–694.
7. A. Tanaka, K. Nitta, R. Maekawa, T. Masuda and T. Higashimura, Polym. J., 1992, 24, 1173–1180.
8. K. Nagai, B. D. Freeman and A. J. Hill, J. Polym. Sci., Part B: Polym. Phys., 2000, 38, 1222–1239.
9. K. D. Dorkenoo and P. H. Pfromm, Macromolecules, 2000, 33, 3747–3751.
10. M. G. Baschetti, M. Ghisellini, M. Quinzi, F. Doghieri, P. Stagnaro, G. Costa and G. C. Sarti, J. Mol. Struct., 2005, 739, 75–86.
11. Y. P. Yampol'Skii, S. M. Shishatskii, V. P. Shantorovich, E. M. Antipov, N. N. Kuzmin, S. V. Rykov, V. L. Khodjaeva and N. A. Platé, J. Appl. Polym. Sci., 1993, 48, 1935–1944.
12. K. Nagai and T. Nakagawa, J. Appl. Polym. Sci., 1994, 54, 1651–1658.
13. L. Starannikova, V. Khodzhaeva and Y. Yampolskii, J. Membr. Sci., 2004, 244, 183–191.
14. K. Takada, H. Matsuya, T. Masuda and T. Higashimura, J. Appl. Polym. Sci., 1985, 30, 1605–1616.
15. T. Masuda, T. Takahashi and T. Higashimura, Macromolecules, 1985, 18, 311–317.
16. Y. Ichiraku, S. A. Stern and T. Nakagawa, J. Membr. Sci., 1987, 34, 5–18.
17. K. Nagai, T. Masuda, T. Nakagawa, B. D. Freeman and I. Pinnau, Prog. Polym. Sci., 2001, 26, 721–798.
18. A. C. Savoca, A. D. Surnamer and C. F. Tien, Macromolecules, 1993, 26, 6211–6216.
19. T. C. Merkel, Z. He, I. Pinnau, B. D. Freeman, P. Meakin and A. J. Hill, Macromolecules, 2003, 36, 6844–6855.
20. S. Matteucci, V. A. Kusuma, D. Sanders, S. Swinnea and B. D. Freeman, J. Membr. Sci., 2008, 307, 196–217.
21. C. H. Lau, P. T. Nguyen, M. R. Hill, A. W. Thornton, K. Konstas, C. M. Doherty, R. J. Mulder, L. Bourgeois, A. C. Y. Liu, D. J. Sprouster, J. P. Sullivan, T. J. Bastow, A. J. Hill, D. L. Gin and R. D. Noble, Angew. Chem., Int. Ed., 2014, 53, 5322–5326.
22. L. Olivieri, S. Ligi, M. G. De Angelis, G. Cucca and A. Pettinau, Ind. Eng. Chem. Res., 2015, 54, 11199–11211.
23. D. Bakhtin, S. Bazhenov, V. Polevaya, E. Grushevenko, S. Makaev, G. Karpacheva, V. Volkov and A. Volkov, Membranes, 2020, 10, 419.
24. G. Golubev, D. Bakhtin, S. Makaev, I. Borisov and A. Volkov, Polymers, 2021, 13, 1922.
25. Q. Fang, Q. Liu, Z. Xie, M. R. Hill and K. Zhang, J. Membr. Sci., 2023, 666, 121139.
26. L. C. Condes, M. T. Webb, T. T. B. Le, W. J. Box, C. M. Doherty, A. Gali, L. Garrido, J. Deng, L. Matesanz-Niño, A. E. Lozano, C. Alvarez, M. Buongiorno Nardelli, A. Striolo, A. J. Hill and M. Galizia, ACS Appl. Mater. Interfaces, 2024, 16, 53843–53854.
27. J. Jia and G. L. Baker, J. Polym. Sci., Part B: Polym. Phys., 1998, 36, 959–968.
28. C. J. Ruud, J. Jia and G. L. Baker, Macromolecules, 2000, 33, 8184–8191.
29. S. D. Bazhenov, I. L. Borisov, D. S. Bakhtin, A. N. Rybakova, V. S. Khotimskiy, S. P. Molchanov and V. V. Volkov, Green Energy Environ., 2016, 1, 235–245.
30. L. Shao, J. Samseth and M.-B. Hägg, J. Membr. Sci., 2009, 326, 285–292.
31. S. D. Kelman, B. W. Rowe, C. W. Bielawski, S. J. Pas, A. J. Hill, D. R. Paul and B. D. Freeman, J. Membr. Sci., 2008, 320, 123–134.
32. D. S. Bakhtin, L. A. Kulikov, S. A. Legkov, V. S. Khotimskiy, I. S. Levin, I. L. Borisov, A. L. Maksimov, V. V. Volkov, E. A. Karakhanov and A. V. Volkov, J. Membr. Sci., 2018, 554, 211–220.
33. B. Lopez-Iglesias, F. Suárez-García, C. Aguilar-Lugo, A. González Ortega, C. Bartolomé, J. M. Martínez-Ilarduya, J. G. De La Campa, Á. E. Lozano and C. Álvarez, ACS Appl. Mater. Interfaces, 2018, 10, 26195–26205.
34. L. C. Condes, A. E. Landry, M. T. Webb, W. J. Box, A. Gali and M. Galizia, J. Membr. Sci., 2024, 709, 123081.
35. D. N. Theodorou, in Materials Science of Membranes for Gas and Vapor Separation, ed. Y. Yampolskii, I. Pinnau and B. Freeman, Wiley, 1st edn, 2006, pp. 49–94.
36. M. Heuchel, D. Fritsch, P. M. Budd, N. B. McKeown and D. Hofmann, J. Membr. Sci., 2008, 318, 84–99.
37. A. Ozcan, R. Semino, G. Maurin and A. O. Yazaydin, Chem. Mater., 2020, 32, 1288–1296.
38. T. T. B. Le, A. Striolo, S. S. Gautam and D. R. Cole, Langmuir, 2017, 33, 11310–11320.
39. T. Bui, A. Phan, D. R. Cole and A. Striolo, J. Phys. Chem. C, 2017, 121, 15675–15686.
40. A. Phan and A. Striolo, Microporous Mesoporous Mater., 2019, 288, 109559.
41. K. Q. Bui, T. T. Bao Le, G. D. Barbosa, D. V. Papavassiliou, S. Razavi and A. Striolo, J. Phys. Chem. Lett., 2024, 15, 8114–8124.
42. C. Divine-Ayela, F. Perez and A. Striolo, Ind. Eng. Chem. Res., 2023, 62, 19893–19906.
43. S. J. Tao, J. Chem. Phys., 1972, 56, 5499–5510.
44. T. L. Dull, W. E. Frieze, D. W. Gidley, J. N. Sun and A. F. Yee, J. Phys. Chem. B, 2001, 105, 4657–4662.
45. A. Stukowski, JOM, 2014, 66, 399–407.
46. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
47. J. Y. Park and D. R. Paul, J. Membr. Sci., 1997, 125, 23–39.
48. J. R. White, C. R. Chim., 2006, 9, 1396–1408.
49. Y. Liu, A. Tang, J. Tan, X. Zhao, C. Chen, D. Wu, Y. Li, P. He and H. Zhang, Polymers, 2020, 12, 2048.
50. D. Hofmann, M. Heuchel, Yu. Yampolskii, V. Khotimskii and V. Shantarovich, Macromolecules, 2002, 35, 2129–2140.
51. W. J. Box, Z. Huang, R. Guo and M. Galizia, J. Membr. Sci., 2022, 656, 120608.
52. J.-H. Zhou, R.-X. Zhu, J.-M. Zhou and M.-B. Chen, Polymer, 2006, 47, 5206–5212.
53. T. C. Merkel, V. Bondar, K. Nagai and B. D. Freeman, J. Polym. Sci., Part B: Polym. Phys., 2000, 38, 273–296.
54. M. Apostolopoulou, R. Dusterhoft, R. Day, M. Stamatakis, M.-O. Coppens and A. Striolo, Int. J. Coal Geol., 2019, 205, 140–154.
55. M. Apostolopoulou, M. Stamatakis, A. Striolo, R. Dusterhoft, R. Hull and R. Day, Transp. Porous Media, 2021, 136, 495–520.
56. D. Argyris, D. R. Cole and A. Striolo, J. Phys. Chem. C, 2009, 113, 19591–19600.
57. T. A. Ho and A. Striolo, J. Chem. Phys., 2013, 138, 054117.
58. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, UK, 2004.
59. R. T. Cygan, V. N. Romanov and E. M. Myshakin, J. Phys. Chem. C, 2012, 116, 13079–13091.
60. M. Galizia, M. G. De Angelis, M. Messori and G. C. Sarti, Ind. Eng. Chem. Res., 2014, 53, 9243–9255.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ta02906f

---