A molecular dynamics simulation study on hydrocarbon ladder polymer membranes for gas separation

Abstract

To address global environmental challenges and support the transition of energy systems, the study of CO₂ capture and separation is at the forefront of scientific research. Utilizing membranes based on polymers of intrinsic microporosity (PIMs) for CO₂ separation presents a promising approach. However, the mechanisms of CO₂ separation in PIMs are not fully understood. In this study, an isobaric model combined with molecular dynamics (MD) simulation was used to explore the adsorptive and diffusive behaviors of CO₂ and N₂ in PIM membranes. We elucidated the gas separation mechanism by analyzing three critical aspects: microporous structure, adsorptive selectivity, and diffusive selectivity. The findings reveal that PIM membranes exhibit advantageous separation characteristics due to their large Brunauer–Emmett–Teller (BET) surface areas and Pore Limiting Diameters (PLDs) that are more compatible with the size of CO₂ molecules. Additionally, the difference in solvation free energy and diffusion rates between the two gases within the membranes significantly contributes to their selectivity. Specifically, CO₂ diffuses within the membrane primarily through a hopping mechanism supplemented by diffusive motion, whereas N₂ relies mainly on diffusion with less hopping. Since dissolution often takes precedence over diffusion in the separation process, it can sometimes lead to less effective diffusion for gas molecules. Moreover, the simulation results indicate that the diffusion behavior of the CO₂/N₂ mixture in PIM membranes is governed by a solubility-driven separation mechanism. This work provides a theoretical foundation for understanding gas transport and separation mechanisms in PIM membranes.

---

1 Introduction

As industrial activities have intensified, the increase in CO₂ emissions has brought numerous environmental challenges, including a rise in global temperatures and strengthening of the greenhouse effect.¹ Notably, this includes sea level rise² and ocean acidification.³ Capturing and storing carbon dioxide from industrial and power generation processes can significantly reduce greenhouse gas emissions in the atmosphere, thereby helping to limit the rise in global average temperatures. For many countries, the separation and capture of CO₂ represent a crucial technological approach to achieving the greenhouse gas emission reduction targets stipulated by the Paris Agreement.

Polymer membrane separation is regarded as an effective method for CO₂ separation. Among all polymer materials, PIMs stand out due to their unique microporous structures. The rigid and contorted skeletons of PIMs hinder chain packing in the solid state, creating significant free volume. These features endow PIMs with superior separation properties compared to traditional polymers,^4–6 particularly for separating CO₂ from other gases such as N₂, H₂, or CH₄.

Due to their outstanding performance, PIMs have established the upper bounds for nearly all gas separations,^7,8 making them a focal point of research in recent years. Since the development of the original PIM-1 prototype, numerous PIMs and their derivatives have been created.^9,10 However, a significant limitation of membrane-based gas separations is the well-known trade-off between permeability and selectivity, which results in performance upper bounds.¹¹ Despite promising advancements, designing highly selective PIMs remains a challenge. Yong et al.¹² were the first to report on the gas separation performance of CO₂/CH₄ using mixed matrix membranes composed of PIM-1 and Matrimid. Since then, PIM-1 has been modified through various techniques aimed at increasing the pore volume and enhancing gas selectivity. These modifications include cyano group post-treatment (such as conversion to carboxylic groups,^13,14 sulfonamides,¹⁵ tetrazoles,^16,17 or amidoxime¹⁸), cross-linking,^19–21 copolymerization,²² development of mixed matrix membranes,^23–28 and the design of specific monomer units.^29–32 Rodriguez et al.³³ reported that at 2 atm, amine-functionalized PIM-1 (PIM-NH2) showed a significant increase in selectivity for CO₂/N₂ mixtures compared to pure gases. Recent research by Lai et al.³⁴ introduced a novel class of carbon-hydrogen ladder polymers that offer both higher selectivity and greater permeability than other membranes in the separation of various industrially relevant gas mixtures while maintaining desirable mechanical and thermal properties in the corresponding films. Unlike traditional PIM membranes, these new polymers opt for more rigid, less flexible bridged-ring structures over highly flexible spirocyclic ones (Fig. 1a). Such PIM membranes promise a combination of mechanical robustness and superior separation efficiency. Experimental conjectures suggest that different hydrocarbon ladder configurations significantly affect chain arrangement and the distribution and connectivity of free volume. Although gas adsorption isotherms and X-ray scattering experiments are commonly used to characterize the porosity of microporous materials, these techniques struggle to capture interactions and relationships at the molecular or atomic level, which are crucial for overcoming separation diffusion bottlenecks. Therefore, employing molecular simulation techniques to reveal the gas separation mechanism of PIM membranes is essential for designing new types of membranes.

Fig. 1 (a) Bridged-ring structure of the DHP monomer and Me2F monomer; (b) spiro-linkage structure of the PIM-1 monomer.

Molecular simulation is a powerful tool for studying adsorptive materials, offering insights into ideal systems and conditions that are experimentally unattainable. For instance, Fang et al.³⁵ conducted simulations on gas permeation through membranes made of intrinsic microporous polymers (PIM-1 and PIM-7). They discovered that nitrile groups in the polymers’ base structure enhanced affinity and increased the fractional free volume (FFV), outperforming most glassy polymer membranes in adsorption and diffusion. Hart et al.³⁶ integrated sulphur into the monomeric framework, offering an additional site, thereby significantly enhancing CO₂ adsorption isotherms and isosteric heat of adsorption in PIMs. Liu et al.³⁷ used molecular dynamics simulations to assess the separation performance of PILP-1 and PILP-3 for CO₂/CH₄ mixtures under constant pressure gradients. This work highlighted the significance of polymer flexibility and membrane plasticization in gas permeation, emphasizing the critical role of gas sorption in CO₂/CH₄ separation. Xu et al.³⁸ examined organic solvent nanofiltration across six PIMs, including the original PIM-1 and five functionalized variants, finding a strong correlation between solvent permeability and the interaction between membrane and solvent properties. These simulation studies offer valuable insights into separation processes in polymer films with intrinsic micropores under specific chemical conditions.

In our study, inspired by previous research, we employed classical molecular dynamics (MD) simulations to explore the gas separation mechanisms of a novel type of hydrocarbon ladder polymer.³⁴ Herein, two types of PIM membranes (Me2F and DHP) were selected as model candidates. Their preparation process is as follows: ladder dinorbornenes were synthesized from 2,7-dibromofluorene derivatives that reacted with excess norbornadiene in the presence of a palladium catalyst; these intermediates subsequently underwent efficient and clean CANAL polymerization with p-dibromo–p-xylene to form ladder polymers Me2F and DHP.³⁴ The structures of Me2F and DHP are shown in Fig. 1a. Unlike classical PIM-1 materials (Fig. 1b), Me2F and DHP feature rigid bridged-ring structures rather than spiro-linkages, limiting the degrees of freedom in polymer chains and inhibiting polymer chain packing. This structural difference provides more accessible channels for gas permeation.

Based on the simulation results, we directly observed the channels and pore structures formed by the rigid stacking of the polymer. We measured the size and distribution of the pores within the membranes as well as the Brunauer–Emmett–Teller (BET) surface area. The entire separation process was analyzed by dividing it into solubility-driven and diffusion-driven components, and we examined the degree of gas dissolution in different membranes through free energy calculations. The kinetic mechanism of gas molecule permeation was investigated using the van Hove function and density functional theory (DFT) calculations. This work provides a foundation for understanding gas separation in these types of membranes and serves as a reference for designing high-performance PIM membranes.

2 Model construction and simulation details

2.1 Model construction

The single polymer chain, comprising 20 monomers of Me2F and DHP, is terminated with hydrogen atoms. Owing to the distinctive rigidity of the new class of carbon–hydrogen ladder PIMs, initial structures were created within a large simulation box, resulting in a polymer system of low density. The insert-molecules module from the Gromacs package³⁹ was utilized to embed 40 polymer chains into the initial simulation box, with dimensions of x = 100 nm, y = 100 nm, and z = 100 nm, as illustrated in Fig. 2a. At present, there is a lack of experimental data regarding the density of Me2F and DHP membranes. In this study, a compression/relaxation strategy consisting of 21 steps was implemented, following the methodology outlined by Colina et al.⁴⁰ It was observed that employing a 21-step gradual decompression scheme ensures that the final simulation density remains unaffected by the maximum pressure applied. The gradual decompression process, involving a maximum temperature of 600 K and a maximum pressure of 50 000 bar, cyclically anneals the system to produce a realistic polymer configuration. This technique has been successfully employed in compressing several PIM membranes. Adhering to the 21-step MD compression and relaxation protocol, the current simulation system was compressed to the target density, resulting in a final box size of x = y = z = 8.9 nm, as depicted in Fig. 2b.

Fig. 2 (a) The initial low-density simulation box; (b) the simulation box after 21-step MD compression; (c) the gas diffusion model of membranes with a constant concentration gradient.

To construct the gas separation system, the post-annealing and compression states were adopted as the initial configuration for developing the CO₂/N₂ mixed gas separation model, resulting in the final separation system as illustrated in Fig. 2c. The simulation box dimensions were set as x = 8.9 nm, y = 8.9 nm, and z = 100 nm. The simulation system was divided into three distinct regions: the left side representing the packing zone, the middle region as the dissolution and permeation zone, and the right side serving as the separation zone.⁴¹ To maintain the volume of the left packing zone constant and prevent the displacement of the polymer film during simulation, a positional constraint was applied to the layer of atoms with the highest Z-coordinate in the PIM membrane.⁴² In Fig. 2c, the purple atoms in the dissolution and permeation zone depict the atoms constrained throughout the simulation process. The volume of the separation zone was slightly larger than that of the packing zone to ensure a lower gas density in the separation zone. To prevent gas molecules from crossing between the packing and separation zones, identical-sized graphene sheets were placed on the left and right sides along the XY-plane. The carbon atoms in the graphene were assigned the charge of zero during the simulation, and uniform interactions between the graphene sheet, membrane, and gas molecules on both sides were maintained. Additionally, a vacuum region was introduced on either side of the graphene sheets to prevent interference between the two graphene layers.³⁸

In the initial phase of the simulation, the system is equilibrated for 100 ns. During the first 50 ns, 100 molecules each of CO₂ and N₂ are introduced to the packing zone, initiating the dissolution of gas molecules into the membrane driven by a concentration gradient. Subsequently, in the following 50 ns, the gas molecules in the separation zone are removed, and new CO₂ and N₂ molecules are inserted into the packing zone to maintain the original number of gas molecules. Following this setup, the balanced system undergoes a 240 ns NVT production process to investigate the separation and diffusion mechanisms of gas molecules within hydrocarbon ladder polymer membranes. The simulation temperature was set as 308 K for this process. This production phase sheds light on the gas separation process and the behaviours of the polymer membrane in gas separation applications.

2.2 Simulation details

In the previous work by Kyle et al.,⁴³ the generalized AMBER force field (GAFF) was employed to model rigid, amorphous, and glassy polymer materials. Their findings indicated that the GAFF could effectively capture crucial details such as pore topology and adsorption capacity. Building on their research, PIMs are also modelled using the GAFF,⁴⁴ while small gas molecules are described with the transferable potential for phase equilibria-united atom (TraPPE-UA) force field.^45–48 For the nonbonded interactions, the Lorentz–Berthelot combination rules for the cross Lennard–Jones (LJ) parameter were adopted. The LJ parameters and atomic information for all gas atoms are provided in Tables S1–S3 (ESI†).

In the simulation process, CO₂ and N₂ molecules were modelled as linear molecules, with virtual atoms introduced to allocate the mass and charge appropriately across the molecules. This approach ensured that the gas molecules had accurate charge distributions within the molecular model. All molecular dynamics simulations were conducted using GROMACS-2021.6-GPU.⁴⁹ The Nosé–Hoover thermostat was utilized to maintain constant temperature, while the Parrinello–Rahman barostat controlled the pressure.⁵⁰ The long-range electrostatic interactions were calculated using the particle-mesh Ewald (PME) method. A cutoff distance of 1.5 nm was used for the LJ potential. The time step for all simulations was set as 1 fs. The final configurations of different systems were visualized by using VMD (version 1.9.3)⁵¹ software.

As the simulation progressed, it was observed that the number of gas molecules in the separation zone increased, leading to a decreased concentration gradient and a reduction in the rate of permeation. To address this issue, a constant concentration differential was maintained throughout the separation and diffusion processes. At regular intervals, denoted by τ, the quantities of gas molecules in both the feed and separation zones were assessed. Gas molecules in the packing zone were replenished according to the initial CO₂ to N₂ ratio (CO₂ = 100, N₂ = 100, and CO₂ : N₂ = 1 : 1). Simultaneously, gas molecules that had permeated into the separation zone were removed. This procedure ensured that the number of CO₂ and N₂ molecules remained constant at 100 in the packing zone, while the separation zone was reset to zero. The original molecular forces and velocities were conserved during the replenishment and extraction processes. Newly added gas molecules underwent a brief dynamic equilibration to adapt to the system conditions. Liu and co-workers⁴² have emphasized the importance of maintaining a constant pressure differential in MD simulations. Their findings suggested that τ values ranging from 150 to 250 ps resulted in stable adsorption levels of CO₂ by using the PIM membrane, without significant fluctuations. To ensure a constant pressure difference between the packing and separation zone, the number of molecules adsorbed by the membrane across different τ values was verified. The results indicated minimal changes in the number of adsorbed gas molecules under varying τ values, as illustrated in Fig. S3 (ESI†). Consequently, the τ interval for gas addition and removal was set at 200 ps in this study.

2.3 Quantum chemistry calculation

To quantitatively investigate the interactions of gas molecules and further explore the distributions of CO₂ and N₂ in hydrocarbon ladder polymers, three interaction models were constructed, which are CO₂–CO₂, N₂–N₂, and CO₂–N₂, respectively. Then, the geometries of all molecules were optimized using density functional theory (DFT) at the GB3LYP-D3/6-31G** level of theory. All calculations are performed using BDF software on the Device Studio platform.^51–56 The interaction region indicator (IRI)⁵⁷ analysis was performed through the Multiwfn (version 3.8_dev) software package.⁵⁸

3 Results and discussion

3.1 Structure and properties of membranes

Molecular simulations offer a unique advantage by providing the atomistic details of the polymeric structure and pore topology.⁵⁹ The model construction scheme involved a deliberate compression over 21 steps, culminating in an 800 ps NPT simulation that provided the converged density of Me2F and DHP polymers at 835.52 kg m⁻³ and 874.21 kg m⁻³, respectively. The equilibrium densities of Me2F and DHP are depicted in Fig. S1 (ESI†).

PIMs are a class of porous glassy polymers known for their substantial free volume, a key structural feature that contributes to their high gas permeability.^16,60 Understanding the pore parameters of these materials is crucial for comprehending their structure and functionality. To evaluate these parameters, Zeo++ software was employed to analyse the pore size metrics of the polymer membrane.⁶¹ This analysis enabled the calculation of several important properties,⁶² including the fractional free volume (FFV), global cavity diameter (GCD), pore limiting diameter (PLD), largest cavity diameter (LCD), BET surface area, and pore size distribution (PSD),⁶³ among others. Table 1 presents detailed BET and FFV data for the membrane, while Fig. 3a and b depict the structure of the membrane's channels and cavities. The entire channel system was divided into two parts: the channel part and the cavity part. Fig. S2 (ESI†) provides a visual representation of the GCD, PLD, and LCD. This comprehensive characterization provides essential insights into the structural attributes of PIMs, which are critical for optimizing their performance in gas separation applications.

Table 1 GCD, PLD, LCD, FFV and BET surface area of the membranes

Mem	GCD (Å)	PLD (Å)	LCD (Å)	FFV (%)	BET surface area (m^2 g^−1) (our work)	BET surface area (m^2 g^−1) (experiment)^34
DHP	11.83 (±1.74)	4.56 (±0.20)	10.28 (±1.24)	53.64	1230.35 (±31.29)	870
Me2F	13.35 (±2.33)	4.12 (±0.17)	11.53 (±1.19)	55.47	1636.96 (±15.34)	1190

---

Fig. 3 The chain distribution and FFV of (a) DHP and (b) Me2F membranes; (c) the pore-size distribution of Me2F and DHP membranes.

In experiments, hydrocarbon ladder polymer membranes are utilized not only for the separation of CO₂ and N₂ but also for other gases such as CH₄ and H₂. To study the pore characteristics, H₂ (with a kinetic radius of 0.14 nm) was selected as a probe molecule. Upon the removal of polymer membranes from all gas molecules, specific information about pore size and surface area can be obtained. The pore size distribution of the Me2F and DHP membranes is shown in Fig. 3c. The minimum pore size of the membranes is approximately 2–3 Å, with the maximum size around 10 Å, and an average size distribution of about 5–6 Å. The distribution indicates that the predominant pore sizes in the Me2F membrane are slightly larger than those in the DHP membrane, which aligns with the observed density distribution. The kinetic diameters of CO₂ and N₂ are 3.47 Å and 3.58 Å, respectively,⁶⁴ making them comparable to the pore sizes in both PIM membranes. This similarity facilitates effective permeation of gases like CO₂ and N₂. It is generally believed that pores with diameters less than 3 Å will be ineffective for gas diffusion, leading to inaccessible areas within the polymer membrane, whereas pores around 10 Å may lack the necessary screening performance. Thus, PIM membranes are expected to be highly effective in separating CO₂ and N₂. Additionally, as seen in Table 1, aside from the pore limiting diameter (PLD), most metrics indicate that the pore sizes in Me2F are slightly larger than those in DHP, which may correlate with the superior permeability observed in Me2F. This finding is consistent with conclusions derived from experimental data (Table S2, ESI†).³⁴ It should be noted that the calculated surface area is greater than the experimental values, potentially due to differences between simulation calculations and experimental measurements.⁶⁵ Furthermore, Me2F has a larger BET surface area compared to DHP, suggesting that gas molecules are more apt to interact with Me2F.

3.2 Absorption behaviours of the CO₂/N₂ mixture in membranes

To explore the gas transport mechanism in a more intuitive manner, we divided the entire separation process into two steps: adsorption and diffusion. Fig. 4a displays the adsorption quantities of CO₂/N₂ mixtures during the permeation process in two distinct membranes. We defined gas molecules as sorbed in the membrane if they were within 0.54 nm of the membrane atoms.⁴² In Fig. 4a, the adsorption of CO₂ is observed to rapidly increase within the initial 10 ns, which subsequently rises gradually from 20 ns to 50 ns, and ultimately stabilizes after 50 ns. This trend suggests that in the initial phase, PIM membranes adsorb CO₂ rapidly. As the adsorption amount increases, the pores of the polymer membrane were occupied progressively by CO₂ molecules, after which the adsorption rate decreases gradually and approaches stability, eventually fluctuating around a specific value. In contrast, the adsorbed quantity of N₂ molecules was significantly lower throughout the entire adsorption process and nearly a constant during the simulation period. We can find that the adsorption quantity of CO₂ is obviously greater than that of N₂.

Fig. 4 (a) The variations of the adsorption number of CO₂ and N₂ molecules versus simulation time; (b) the variation of the membrane density versus simulation time; solvation free energies of two gas molecules in (c) DHP and (d) Me2F membranes.

This disparity in adsorption quantities indicates that CO₂ has a much stronger affinity for the membrane than N₂. This can be attributed to the molecular interactions and the more fitting kinetic diameter of CO₂ relative to the pore sizes in the PIM membranes, which enhances CO₂ adsorption. This strong preferential adsorption of CO₂ is crucial for effective gas separation, as it ensures selective adsorption followed by diffusion through the membrane, thereby optimizing separation efficiency.

In our analysis of the separation process, we focused primarily on the solvation aspect. Fig. 4b illustrates the localized density distribution, showing that PIM membranes experience some swelling after gas adsorption. This swelling leads to a reduction in density from the initial dense regions within the simulation model, indicating that the membrane structure dynamically adjusts in response to the permeation process. Eventually, the density stabilizes at a constant value, signalling that the adsorption of gas molecules has reached a saturated state.

To further explore the solvated adsorption process, we calculated the solvation free energies for the two gases in different membranes. This involved computing the change in free energy (ΔG) between the state of the membrane containing gas molecules (state A) and the pure membrane state (state B) as a function of the coupling parameter λ, which represents the transition extent from state A to state B. As shown in Fig. 4c and d, in DHP, ΔG_DHP–CO2 is −10.21 kJ mol⁻¹ and ΔG_DHP–N2 is −3.87 kJ mol⁻¹, whereas in Me2F, ΔG_Me2F–CO2 is −9.48 kJ mol⁻¹ and ΔG_Me2F–N2 is −3.58 kJ mol⁻¹. The significant differences in free energy between CO₂ and N₂ for both DHP and Me2F membranes indicate much stronger adsorption of CO₂ compared to N₂. These differences in solvation free energies facilitate the separation process in these membranes. It is well-known that the performance of gas separation membranes is predominantly influenced by two critical parameters: the solubility coefficient (S) of gases within the membrane material and the diffusion coefficient (D) of gas molecules through the membrane matrix. Together, these parameters determine the efficiency and selectivity of the separation process. They are pivotal indicators for evaluating and optimizing the performance of gas separation membranes, guiding both design improvements and operational adjustments.

Adsorption selectivity is a criterion for evaluating adsorption priority in binary mixed systems:⁶⁶

The value of adsorption selectivity S_i/j larger than 1 indicates that component i is preferentially adsorbed. The x_i is the mole fraction of gas i in the adsorbed phase, and y_i is the mole fraction of gas i in the membrane bulk phase. To study the competitive adsorption behaviour of CO₂ and N₂, the adsorption selectivity of CO₂/N₂ is defined as follows:

where S_CO2/N2 is the adsorption selectivity of CO₂ over N₂. The partial pressure may play an important role in the competitive adsorption behaviour.⁶⁷ In our work, because the ratio of CO₂ to N₂ in the packing zone is always 1 : 1, so there is no partial pressure. According to eqn (2), in DHP and Me2F membranes, the solubility selectivity is 5.30 (±0.27) and 4.71 (±0.08), respectively.

3.3 Permeation behaviours of the CO₂/N₂ mixture through membranes

To assess the impact of gas diffusion on the performance of separation membranes, we analysed the dynamic behaviours involved in the diffusion process. By accurately measuring the diffusion rates and interaction conditions of gases within the membrane materials, we can optimize membrane structures and improve their practical performance.

Fig. 5a illustrates the relationship between the permeation flux and diffusion time. Initially, N₂ diffusion efficiency is higher than that of CO₂ because the voids in the PIMs are not yet occupied by CO₂. In this early phase, CO₂ molecules are adsorbed within the membrane cavities and cannot diffuse readily, whereas N₂, experiencing weaker interactions with the membrane, moves more freely. As the simulation progresses, CO₂ adsorption reaches saturation (see Fig. 4a), fully occupying the membrane's voids. After further gas adsorption, the permeability efficiency of CO₂ was significantly higher than that of N₂, showing a linear increase. Furthermore, the Me2F membrane demonstrates a superior permeation rate compared to DHP. The rate of gas migration serves as an effective metric for gauging the diffusion velocity of gas molecules moving within a membrane. This parameter reflects the dynamic characteristics of the movement of gas molecules within the separation membrane medium, offering a quantitative basis for evaluating the efficiency of gas diffusion through the membrane. The migration rate of gas molecules in membranes can be determined by the mean square displacement function (MSD), which can be expressed as:

where N represents the total number of gas molecules, and r_k(t) denotes the position of the molecule k at time t. MSD calculations were based on the ensemble average of the molecular trajectories using eqn (3), and the results are presented in Fig. 5b. The Einstein equation is employed to derive the diffusion coefficients:

Fig. 5 (a) Permeation numbers of CO₂ and N₂ molecules versus simulation time in Me2F and DHP membranes. (b) Mean square displacement (MSD) of gas molecules in membranes. (c) The van Hove function for CO₂ in DHP; (d) N₂ in DHP; (e) CO₂ in Me2F; (f) N₂ in Me2F.

MDAnalysis⁶⁸ was used to analyse the MSD of gas molecules within the system, allowing us to calculate their diffusion coefficients. Table 2 presents these diffusion coefficients for gas molecules at 308 K, as determined by eqn (4). These values are averaged from three independent simulations. Fig. 5b shows that the permeability coefficient of CO₂ is significantly lower than that of N₂, indicating a much lower migration rate for CO₂ within the membrane. Based on this, one might consider that the membrane's selectivity favours N₂ over CO₂. However, both experimental and simulation results reveal the opposite: the membrane demonstrates higher selectivity for CO₂.³⁴

Table 2 The diffusion coefficients of CO₂ and N₂ in DHP and Me2F membranes at 308 K

Gas	DHP (10^−5 cm^2 s^−1)	Me2F (10^−5 cm^2 s^−1)
CO2	301.44 (±9.13)	424.77 (±60.95)
N2	703.99 (±67.68)	758.66 (±41.25)

---

To understand the unexpected appearance and further explore the transport mechanisms of gas molecules within the membrane, their transport behaviours were analysed using the van Hove correlation function, as shown in Fig. 5c–f. This analysis revealed that CO₂ and N₂ molecules exhibit distinct diffusion mechanisms in the two types of membranes. The van Hove function displays two distinct peaks for both gases. The first peak corresponds to the local motion of gas molecules, while the second peak indicates their diffusion motion over longer distances. For CO₂, the first peak is located at 0.1 nm, suggesting that CO₂'s movement is significantly restricted within the membrane structure. In contrast, while N₂ also shows a characteristic peak at 0.1 nm, its intensity is much lower than that of CO₂, indicating that N₂'s local motion is less constrained. The second peak, associated with long-range movements, is located at greater distances and has broader shape. As simulation time progresses, these peaks shift further outward, illustrating the transition of gas molecules from localized movements to long-range diffusion. Comparing the second peaks of CO₂ and N₂ reveals that N₂ experiences more dominant long-range movements, which aligns well with the observed diffusion results. This suggests that despite CO₂'s higher selectivity in the membrane, N₂ has a greater propensity for long-range diffusion, which is in good agreement with our diffusion results.

Comparing the van Hove correlation function for the same gas species across different types of separation membranes reveals notable differences in their transport behaviors. For CO₂, the first peak in DHP membranes is more pronounced compared to Me2F, indicating that CO₂ experiences greater constraints within DHP. This suggests that because of these limitations, CO₂ in DHP is more difficult to diffuse out of the membrane than CO₂ in Me2F. Similarly, the diffusion behavior of N₂ across both membrane types follows a pattern akin to CO₂. DHP imposes more significant restrictions on N₂, whereas in Me2F, N₂ experiences more extensive jumps. This results in a more pronounced diffusion effect for N₂ in Me2F compared to DHP.

Within the same type of membrane, N₂ experiences fewer restrictions compared to CO₂ and primarily follows classical Fickian diffusion kinetics. This suggests that N₂'s diffusion behavior is dominated by free diffusion, with hopping mechanisms playing a secondary role. Conversely, CO₂ moves shorter distances over the same period due to stronger constraints, indicating that for CO₂, the hopping mechanism is primary, while diffusion is secondary. This highlights that the two gases exhibit distinct transport mechanisms. In both types of membranes, Me2F demonstrates higher permeation flux as it imposes relatively weaker constraints on gas molecules, as shown in Fig. 5a.

To gain insights into the reasons behind these different diffusion mechanisms, the RDF was used to analyze the interactions between gas molecules and the PIMs, as well as the changes in coordination numbers of the gas molecules dispersed within the membranes. Fig. 6 shows that the RDFs of these two gas molecules in the two membranes have the same shape. Compared to Me2F, DHP exhibits the peak at the same position, indicating that both membranes have a close degree of coordination tightness for gases. By integrating the area under the RDF curve within the shell radius, the coordination number of CO₂ and N₂ molecules in the DHP membrane is 2.47 and 1.17, and 2.65 and 1.73 in the Me2F membrane, respectively. The coordination numbers of the two gases are obviously different in different membranes, which reflects that the both membranes can better separate CO₂/N₂.

Fig. 6 RDFs (solid lines) and CNs (dashed lines) between the characteristic atom CX and gas molecules in (a) DHP and (b) Me2F membranes.

To further compare the differences between the two models, we calculated the residence times of gas molecules. The auto-correlation function (ACF) of the residence time can represent the probability of a group of particles staying in the specified region or the survival probability (SP). The decay speed can reflect the kinetic state of the particles by eqn (5), where t is the time step, N(t₀) is the number of gas molecules in the specified region at time t₀, and N(t₀,t₀ + t) is the number of gas molecules in each frame between t₀ and t₀ + t. The angle bracket indicates that the ensemble average is taken for the initial value t₀ at all times.

Using eqn (6) to fit the ACF, the average residence time of the particles in the specified region can be obtained.

where τ represents the retention time, and β denotes the fitting coefficient.

Fig. 7 illustrates the survival probability of gas molecules within the membrane. In DHP, the average retention time for CO₂ and N₂ is 18.10 ps and 11.17 ps, respectively. In contrast, the retention time increases to 26.52 ps for CO₂ and 13.89 ps for N₂ in Me2F. These data indicate that CO₂'s movement is more restricted compared to N₂, making CO₂ permeation through these membranes more challenging.

Fig. 7 Residence autocorrelation function of the retention time of gas molecules in the PIM membrane.

In order to quantify the interaction between gas molecules and membranes, the interaction energy is calculated and shown in Fig. 8. We can see that the van der Waals and coulomb interactions between CO₂ and the membrane are much larger than that of N₂, and the membrane is more selective to CO₂. This finding is consistent with the experimental result.³⁴ Based on the above simulation results of the gas adsorption process, we can obtain that the adsorption effect of CO₂ is stronger than that of N₂. That is say, the diffusion effect was weaker than that of N₂. But a contrary conclusion between experimental and simulation results was obtained, where the permeation effect of CO₂ is greater than that of N₂. To unveil this contradictory phenomenon, we calculated the binding strength between gas molecules by quantum chemistry calculations. As shown in Fig. 9, the binding energies of CO₂–CO₂, N₂–N₂ and CO₂–N₂ are −2.47, −0.71, and −1.40 kcal mol⁻¹, respectively. The interaction region indicator (IRI) analysis (Fig. 10) shows that the binding energy is mainly derived from the van der Waals interaction between carbon and oxygen atoms. According to the binding strength order, we can infer that CO₂ molecules tend to form CO₂ clusters, and CO₂ molecules are not transported as a single molecule in the membrane during diffusion separation. Through the trajectory movie shown in Fig. S4 (ESI†), we preliminarily confirm the existence of CO₂ molecular clusters. It also reflects the reason why CO₂ is more restricted in the membrane.

Fig. 8 van der Waals interaction and coulomb interaction energy between gas and membranes.

Fig. 9 Binding energy of CO₂–CO₂, N₂–N₂ and CO₂–N₂.

Fig. 10 The interaction region indicator analysis of the CO₂ dimer.

To address the discrepancy between the observed diffusion rates and selectivity of gases, it's important to consider both adsorption and diffusion kinetics. Although N₂ diffuses faster, CO₂ exhibits greater selectivity due to its stronger interaction with the membrane, which affects both adsorption and diffusion processes. Based on the solution–diffusion mechanism, the permeability P in polymer membranes can be expressed as:

P = S × D (7)

The separation factor between two substances i and j in the membrane is typically assessed through the ideal permeation selectivity.³⁵

a_i/j = P_i/j = S_i/j × D_i/j (8)

where P_i/j is the selectivity for both i and j gases, S_i/j is the solubility selectivity, and D_i/j is the diffusivity selectivity. Table 3 indicates that CO₂ has stronger solubility selectivity over N₂ in both membranes, and the diffusion selectivity is lower than N₂. This permeation mechanism elucidated the interesting simulated phenomenon mentioned above, where the gas with faster diffusion paradoxically exhibits slower transport. It can be concluded that in these separation membranes, the separation efficiency is codetermined by solubility and diffusion. Even though the diffusion process favors N₂, the overall selectivity is weakened for N₂ due to competition between dissolution and diffusion factors. The selectivity obtained by our work is similar to the experimental conclusion, that is, the selectivity of CO₂ is much higher than that of N₂.

Table 3 Comparison of gas permeation properties of polymers (DHP and Me2F)

i = CO2, j = N2	S i/Sj	D i/Dj	P i/Pj
DHP	5.30 (±0.27)	0.43 (±0.01)	2.27 (±0.18)
Me2F	4.71 (±0.08)	0.56 (±0.05)	2.63(±0.29)

---

4 Conclusion and outlook

In this work, molecular dynamics simulations were conducted to investigate the separation of CO₂/N₂ mixtures in two types of microporous hydrocarbon ladder polymer membranes. The mechanism of gas separation was revealed by analysing the pore size, structure and kinetic and thermodynamic data of the polymers. According to the separation process, we divide the whole process into two key stages: adsorption and separation. During the initial adsorption stage, CO₂ gas rapidly diffused into the membrane and attained equilibrium, while N₂ maintained a relatively lower permeation level. Once CO₂ permeation reached equilibrium, its diffusion efficiency will exceed that of N₂. Both membranes demonstrated high selectivity for CO₂ gases, primarily attributed to solubility selectivity. In the second stage, CO₂ predominantly utilized a hopping mechanism and secondary diffusion mechanism to permeate through the membrane. However, N₂ mainly followed a traditional diffusion mechanism and an auxiliary hopping mechanism. Moreover, through the calculation of interaction energy between gas molecules and membranes, we found that CO₂ tends to form clusters in membranes. Finally, the selectivity of Me2F to CO₂/N₂ was 2.63, and that of DHP to CO₂/N₂ was 2.27. Both membranes showed greater solubility selectivity than diffusion selectivity. That is, CO₂/N₂ through new PIMs mainly conforms to the solution–selective separation mechanism. The current work provides a theoretical basis for the design and manufacture of high-performance membranes.

In this work, the absolute values of the ideal selectivity in the experiment and simulation have a certain deviation (Table S4 and Fig. S5, ESI†). The discrepancy may primarily stem from variations in the driving pressure applied in experiments and simulations. Additionally, potential inaccuracies in the force field parameters cannot be entirely ruled out. In experiments, PIM membranes are inherently heterogeneous, whereas our simulations utilized periodic boundary conditions, simplifying the intricate and multiscale structure of real membranes. Despite efforts to replicate experimental conditions during model construction, the idealized treatment of the PIM membrane model likely contributed to these inconsistencies. The accuracy of the force field is paramount for accurately describing intermolecular interactions, particularly when applied to heterogeneous materials like PIM membranes. Future research could benefit from exploring machine learning-based force fields (MLFFs) specifically tailored for PIM membranes, alongside advanced modeling techniques, to enhance simulation accuracy. MLFFs, trained on high-precision quantum chemical data, offer the potential to better describe weak interactions and adapt to the complexity of heterogeneous materials. Implementing these advancements could significantly improve the precision and reliability of PIM membrane simulations, paving the way for more accurate predictive models in membrane science.

Author contributions

W. X. Tian: data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft, writing – review & editing. L. D. Gong: conceptualization, data curation, formal analysis, investigation, methodology, software, project administration, resources, supervision, writing – review & editing. C. Y. Yu: conceptualization, data curation, formal analysis, investigation, methodology, software, project administration, resources, supervision, writing – review & editing. Y. F. Zhou: conceptualization, project administration, resources, supervision, writing – review & editing.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge the financial support from the National Natural Science Foundation of China (No. 22377048). We gratefully acknowledge HZWTECH for providing computation facilities.

Notes and references

1. L. Al-Ghussain, Environ. Prog. Sustainable Energy, 2018, 38, 13–21.
2. A. Dutton, A. E. Carlson, A. J. Long, G. A. Milne, P. U. Clark, R. DeConto, B. P. Horton, S. Rahmstorf and M. E. Raymo, Science, 2015, 349, aaa4019.
3. K. J. Kroeker, R. L. Kordas, R. Crim, I. E. Hendriks, L. Ramajo, G. S. Singh, C. M. Duarte and J. P. Gattuso, Global Change Biol., 2013, 19, 1884–1896.
4. C. G. Bezzu, M. Carta, A. Tonkins, J. C. Jansen, P. Bernardo, F. Bazzarelli and N. B. McKeown, Adv. Mater., 2012, 24, 5930–5933.
5. M. Carta, R. Malpass-Evans, M. Croad, Y. Rogan, J. C. Jansen, P. Bernardo, F. Bazzarelli and N. B. McKeown, Science, 2013, 339, 303–307.
6. B. S. Ghanem, R. Swaidan, X. Ma, E. Litwiller and I. Pinnau, Adv. Mater., 2014, 26, 6696–6700.
7. B. Comesaña-Gándara, J. Chen, C. G. Bezzu, M. Carta, I. Rose, M.-C. Ferrari, E. Esposito, A. Fuoco, J. C. Jansen and N. B. McKeown, Energy Environ. Sci., 2019, 12, 2733–2740.
8. R. Swaidan, B. Ghanem and I. Pinnau, ACS Macro Lett., 2015, 4, 947–951.
9. P. M. Budd, B. S. Ghanem, S. Makhseed, N. B. McKeown, K. J. Msayib and C. E. Tattershall, Chem. Commun., 2004, 230–231.
10. P. M. Budd, E. S. Elabas, B. S. Ghanem, S. Makhseed, N. B. McKeown, K. J. Msayib, C. E. Tattershall and D. Wang, Adv. Mater., 2004, 16, 456–459.
11. L. M. Robeson, J. Membr. Sci., 2008, 320, 390–400.
12. W. F. Yong, F. Y. Li, Y. C. Xiao, P. Li, K. P. Pramoda, Y. W. Tong and T. S. Chung, J. Membr. Sci., 2012, 407–408, 47–57.
13. N. Du, G. P. Robertson, J. Song, I. Pinnau and M. D. Guiver, Macromolecules, 2009, 42, 6038–6043.
14. J. Weber, N. Du and M. D. Guiver, Macromolecules, 2011, 44, 1763–1767.
15. C. R. Mason, L. Maynard-Atem, N. M. Al-Harbi, P. M. Budd, P. Bernardo, F. Bazzarelli, G. Clarizia and J. C. Jansen, Macromolecules, 2011, 44, 6471–6479.
16. N. Du, H. B. Park, G. P. Robertson, M. M. Dal-Cin, T. Visser, L. Scoles and M. D. Guiver, Nat. Mater., 2011, 10, 372–375.
17. N. Du, G. P. Robertson, M. M. Dal-Cin, L. Scoles and M. D. Guiver, Polymer, 2012, 53, 4367–4372.
18. H. A. Patel and C. T. Yavuz, Chem. Commun., 2012, 48, 9989–9991.
19. N. Du, M. M. D. Cin, I. Pinnau, A. Nicalek, G. P. Robertson and M. D. Guiver, Macromol. Rapid Commun., 2011, 32, 631–636.
20. N. Du, M. M. Dal-Cin, G. P. Robertson and M. D. Guiver, Macromolecules, 2012, 45, 5134–5139.
21. F. Y. Li, Y. Xiao, T.-S. Chung and S. Kawi, Macromolecules, 2012, 45, 1427–1437.
22. T. Emmler, K. Heinrich, D. Fritsch, P. M. Budd, N. Chaukura, D. Ehlers, K. Rätzke and F. Faupel, Macromolecules, 2010, 43, 6075–6084.
23. J. Ahn, W.-J. Chung, I. Pinnau, J. Song, N. Du, G. P. Robertson and M. D. Guiver, J. Membr. Sci., 2010, 346, 280–287.
24. A. F. Bushell, M. P. Attfield, C. R. Mason, P. M. Budd, Y. Yampolskii, L. Starannikova, A. Rebrov, F. Bazzarelli, P. Bernardo, J. Carolus Jansen, M. Lanč, K. Friess, V. Shantarovich, V. Gustov and V. Isaeva, J. Membr. Sci., 2013, 427, 48–62.
25. A. F. Bushell, P. M. Budd, M. P. Attfield, J. T. A. Jones, T. Hasell, A. I. Cooper, P. Bernardo, F. Bazzarelli, G. Clarizia and J. C. Jansen, Angew. Chem., Int. Ed., 2012, 52, 1253–1256.
26. C. Das Neves Gomes, O. Jacquet, C. Villiers, P. Thuéry, M. Ephritikhine and T. Cantat, Angew. Chem., Int. Ed., 2011, 51, 1.
27. M. G. De Angelis, R. Gaddoni and G. C. Sarti, Ind. Eng. Chem. Res., 2013, 52, 10506–10520.
28. M. M. Khan, V. Filiz, G. Bengtson, S. Shishatskiy, M. M. Rahman, J. Lillepaerg and V. Abetz, J. Membr. Sci., 2013, 436, 109–120.
29. M. Carta, K. J. Msayib, P. M. Budd and N. B. McKeown, Org. Lett., 2008, 10, 2641–2643.
30. D. Fritsch, G. Bengtson, M. Carta and N. B. McKeown, Macromol. Chem. Phys., 2011, 212, 1137–1146.
31. B. S. Ghanem, Polym. Chem., 2012, 3, 96–98.
32. S. Li, H. J. Jo, S. H. Han, C. H. Park, S. Kim, P. M. Budd and Y. M. Lee, J. Membr. Sci., 2013, 434, 137–147.
33. K. Mizrahi Rodriguez, F. M. Benedetti, N. Roy, A. X. Wu and Z. P. Smith, J. Mater. Chem. A, 2021, 9, 23631–23642.
34. H. W. H. Lai, F. M. Benedetti, J. M. Ahn, A. M. Robinson, Y. Wang, I. Pinnau, Z. P. Smith and Y. Xia, Science, 2022, 375, 1390–1392.
35. W. Fang, L. Zhang and J. Jiang, Mol. Simul., 2010, 36, 992–1003.
36. K. E. Hart, L. J. Abbott, N. B. McKeown and C. M. Colina, Macromolecules, 2013, 46, 5371–5380.
37. J. Liu and J. Jiang, J. Phys. Chem. C, 2019, 123, 6607–6615.
38. Z. Zhao and J. Jiang, J. Membr. Sci., 2020, 608, 118173.
39. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
40. G. S. Larsen, P. Lin, K. E. Hart and C. M. Colina, Macromolecules, 2011, 44, 6944–6951.
41. B. Ma, C. Chen, Z.-X. Du and L. Zhang, J. Mol. Liq., 2023, 386, 122362.
42. X. Kong and J. Liu, J. Phys. Chem. C, 2019, 123, 15113–15121.
43. K. E. Hart, L. J. Abbott and C. M. Colina, Mol. Simul., 2013, 39, 397–404.
44. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
45. N. Lubna, G. Kamath, J. J. Potoff, N. Rai and J. I. Siepmann, J. Phys. Chem. B, 2005, 109, 24100–24107.
46. M. G. Martin and J. I. Siepmann, J. Phys. Chem. B, 1998, 102, 2569–2577.
47. N. Rai and J. I. Siepmann, J. Phys. Chem. B, 2007, 111, 10790–10799.
48. C. D. Wick, J. M. Stubbs, N. Rai and J. I. Siepmann, J. Phys. Chem. B, 2005, 109, 18974–18982.
49. S. Nosé, Mol. Phys., 2006, 52, 255–268.
50. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
51. W. Liu, G. Hong, D. Dai, L. Li and M. Dolg, Theor. Chem. Acc., 1997, 96, 75–83.
52. W. Liu, F. Wang and L. Li, Recent Advances in Relativistic Molecular Theory, World Scientific, Singapore, 2004, vol. 5, pp. 257–282.
53. W. Liu, F. Wang and L. Li, J. Theor. Comput. Chem., 2011, 02, 257–272.
54. Z. Wang, Z. Li, Y. Zhang and W. Liu, J. Chem. Phys., 2020, 153, 164109.
55. Y. Zhang, B. Suo, Z. Wang, N. Zhang, Z. Li, Y. Lei, W. Zou, J. Gao, D. Peng, Z. Pu, Y. Xiao, Q. Sun, F. Wang, Y. Ma, X. Wang, Y. Guo and W. Liu, J. Chem. Phys., 2020, 152, 064113.
56. Hongzhiwei Technology, Device Studio, Version 2023A, China, 2023, Available online: https://cloud.hzwtech.com/web/pro-ductservice?id=6 (accessed on Aug., 23rd).
57. T. Lu and Q. Chen, Chem.: Methods, 2021, 1, 231–239.
58. T. Lu and F. Chen, J. Comput. Chem., 2011, 33, 580–592.
59. K. E. Jelfs and A. I. Cooper, Curr. Opin. Solid State Mater. Sci., 2013, 17, 19–30.
60. N. B. McKeown and P. M. Budd, Macromolecules, 2010, 43, 5163–5176.
61. T. F. Willems, C. H. Rycroft, M. Kazi, J. C. Meza and M. Haranczyk, Microporous Mesoporous Mater., 2012, 149, 134–141.
62. S. Matteucci, Y. Yampolskii, B. D. Freeman and I. Pinnau, Mater. Sci. Membr. Gas Vap. Sep., 2006, 1–47.
63. L. Sarkisov and A. Harrison, Mol. Simul., 2011, 37, 1248–1257.
64. N. Mehio, S. Dai and D.-E. Jiang, J. Phys. Chem. A, 2014, 118, 1150–1154.
65. D. A. Gómez-Gualdrón, P. Z. Moghadam, J. T. Hupp, O. K. Farha and R. Q. Snurr, J. Am. Chem. Soc., 2016, 138, 215–224.
66. L. Brochard, M. Vandamme, R. J. Pellenq and T. Fen-Chong, Langmuir, 2012, 28, 2659–2670.
67. J. Sun, C. Chen, Y. Zhang, W. Li and Y. Song, Chem. Eng. J., 2022, 430, 133172.
68. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.

---

Footnote

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

---