MOF surface morphology governs interfacial pore architecture and CO₂ dynamics in mixed matrix membranes

Abstract

Mixed matrix membranes (MMMs), which embed metal–organic frameworks (MOFs) within polymers, offer a promising platform for next-generation, energy-efficient separations. However, the nano-structuring of the MOF/polymer interface and its influence on the MMM performance remains poorly understood. Here, we uncover two fundamental design principles that bridge this gap enabled by an automated, graph theory enhanced molecular simulation platform. First, we demonstrate that MOF surface morphology, specifically its planarity and roughness, plays a decisive role in shaping the topology of the interfacial pore network, including its dimensionality, connectivity, and spatial organization. Second, we show that this pore topology critically governs interfacial CO₂ dynamics: highly interconnected and continuous networks facilitate efficient translational and rotational motion, whereas fragmented architectures severely limit molecular mobility. Beyond providing a deep molecular-level understanding, this work introduces a new design paradigm: deliberate tuning of MOF surface morphology emerges as a powerful strategy to control interfacial nanostructure and optimize gas dynamics. Together, these findings open an unexplored pathway for the rational design of high-performance MMMs for advancing energy-efficient separation technologies.

---

1. Introduction

The global imperative to reduce energy consumption and mitigate environmental impacts has spurred a concerted effort to develop more energy-efficient and environmentally sustainable separation technologies. Traditional methods such as distillation are known for their high energy requirements, prompting a critical need for greener alternatives.¹ Membrane-based separation has emerged as an effective technology in this endeavor, offering low-energy solutions for a wide range of liquid and gas mixtures.² However, most of the first-generation polymeric membranes face a trade-off between permeability and selectivity.³ Typically, more permeable polymers tend to be less selective, and vice versa. This phenomenon was evidenced by the Robeson upper bound.^4,5 To overcome this challenge, a new class of hybrid membranes known as Mixed Matrix Membranes (MMMs) has garnered significant attention.^6,7 By incorporating highly selective ordered porous materials as fillers i.e. zeolites, metal–organic frameworks (MOFs), covalent–organic frameworks (COFs) and carbon nanotubes, into polymeric matrices, MMMs combine the advantages of both polymers (mechanical stability, facile processing, low cost) and inorganic/organic/hybrid materials (high selectivity), presenting a promising avenue for various gas separation applications.^8–11 This holds particularly true for MMMs combining MOFs with a wide set of polymers that have shown attractive separation performance for a range of gas mixtures.¹²

These composite systems imply the creation of a boundary zone between the polymer matrix and the dispersed filler. This MOF/polymer interface has been shown to be pivotal in ensuring the fabrication of effective MMMs for molecular separation.^6,13–24 Precise control of interfacial properties is critical: (1) strong MOF/polymer adhesion prevents defects formation and preserves membrane stability;²⁵ (2) uniform MOF dispersion within the polymer matrix avoids aggregation and guest transport barriers;²⁶ and (3) an optimized interfacial pore network favors efficient guest molecule transfer from the polymer matrix into MOF selective pores.^21,22 Therefore, a systematic study is required to deepen our understanding of the MOF/polymer interfacial pore structuring and its impact on MMMs performance to further identify emergent design rules of optimal MMMs. The interfacial pore structuring has not yet been sufficiently investigated, largely due to the complexity of its characterization from both experimental and modelling perspectives.

Our group has deployed a computational strategy integrating force field equilibrium Molecular Dynamics (FF-EMD) and quantum mechanical (QM) calculations¹³ to deliver first insight into this fundamental question at the microscopic level with the exploration of a series of MOF/polymer composites.^{16,21,22,27–29} We demonstrated that a fine-tuning of the interfacial MOF/polymer interactions by considering flexible and/or functionalized polymers,^14,17 or editing structural defects at the MOF surface^16,30 govern the adhesion between the 2 components to ensure the processability of continuous and mechanically stable membranes. Our most recent study on MMMs constructed by assembling the square-shaped ultra-small pore AlFFIVE-1-Ni (KAUST-8) and two distinct polymers, i.e. PIM-1 and 6FDA-DAM revealed that the interfacial pore shape plays also a key role for an optimal molecular transport in the overall MMM.^21,22 The predicted diffuse interfacial pore shape of the AlFFIVE-1-Ni/PIM-1 composite was demonstrated to enable a more effective CO₂ migration from the polymer region towards the entrance of the selective MOF filler. This prediction was further supported by the excellent CO₂ permeability shown by the corresponding [001]-oriented nanosheets AlFFIVE-1-Ni/PIM-1 MMM. The selection of this prototypical ultra-small pore MOF delivered a compelling proof of concept, demonstrating the potential of modulating MOF/polymer interfacial pore structuring to enhance the molecular transport throughout the overall highly orientated MMM and hence its gas permeability performance.

To generalize this overall concept, a systematic assessment of the role played by distinct nature of interfacial pore structuring (dimension, shape, interconnectivity) on the MOF/polymer adhesion and interfacial molecular diffusion is required. On this path, we selected a set of CO₂ selective narrow-pore MOFs sorbents in addition to the square-channel-like AlFFIVE-1-Ni. This includes CALF-20,²⁸ MIL-53-NH₂ (ref. 31 and 32) and Zr-fum-fcu,²⁹ which feature distinct pore aperture shapes, oval, lozenge, and triangular respectively. We specifically chose exposed surfaces oriented perpendicular to the pore apertures and spanning a broad range of surface planarity/roughness to interrogate how morphology governs interfacial pore-network topology, i.e. [001] (AlFFIVE-1-Ni), [100] CALF-20, [001] MIL-53-NH₂ and [110] Zr-fum-fcu facets. The Delaunay tessellation approach as implemented in Qhull³³ was applied to precisely analyze the roughness and morphologies of these four distinct MOF surfaces by assessing their respective planarity deviation (P_dev) (P_dev ≈ 1: highest planarity, lowest P_dev indicates increased surface irregularity). Fig. 1 highlights the deliberate selection of these four MOF surfaces spanning a range of planarity variations, providing a systematic framework to investigate how surface roughness governs the interfacial pore network formation upon integration with a polymer.

Fig. 1 Illustration of the exposed MOF slab models orientated perpendicularly to the MOF pore entrance: [001] ALFFIVE-1-Ni, [100] CALF-20, [001] MIL53-NH₂ and [110] Zr-fum-fcu. (A) Top view and (B) side view of the respective MOFs including Delaunay tessellation surface. (C) 3D Delaunay tessellation plane representation of the MOF surface and the corresponding planarity deviation metrics. Herein, the 3D Voronoi diagram was constructed based on the MOF surface atoms with a grid of 1 Å considered to detect atoms closest to the surface within a 2 Å radius. Due to the particularly large pore aperture of Zr-fum-fcu, this radius was extended up to 4 Å.

Composites were constructed by integrating these different MOF surfaces with 6FDA-DAM, a highly permeable perfluorinated dianhydride-based polyimide, chosen as a model polymer to generate distinct nature of interfacial pore structuring. Notably, this polymer is widely used in MMMs for its high thermal and chemical stability, good mechanical strength, and excellent processability.^21,22 To that purpose we deployed our well-trained computational approach,¹³ based on a combination of FF-EMD and QM calculations. Notably, our previous workflow analyses provided local or ensemble-average descriptors of the interface, but it was not possible to characterize the three-dimensional connectivity of the interfacial free volume: whether voids are isolated, weakly connected, or form continuous multi-path channels. Because interfacial transport depends on percolating pathways and bottlenecks rather than only on average free volume, a topology-aware description is required. Herein, we established an automated graph theory enhanced framework that expands our computational capabilities by making high-resolution analysis of complex interfacial pore networks possible as schematized in Fig. 2. This decisive implementation especially enabled to convert the accessible interfacial free volume (defined by our Delaunay-based segmentation) into a weighted network whose nodes are accessible grid points and whose edges encode spatial proximity. From this representation we can extract three complementary descriptors: eccentricity (spatial reach/longest path within a pore domain), betweenness centrality (bottlenecks and articulation points that control through-interface passage), and assortativity (tendency of similar void regions to connect, i.e., cohesiveness domain). These metrics are insensitive to arbitrary 2D slicing and capture full 3D connectivity, thereby complementing conventional density, RDF and PSD analyses. This overall upgraded platform enables a systematic analysis of the interfacial pore structuring within the constructed MOF/6FDA-DAM composites with unparalleled precision. We demonstrate that the exposed MOF surface morphology governs the interfacial pore architecture, with variations in MOF surface roughness leading to the formation of distinct pore networks exhibiting different structural and textural characteristics. These interfacial features, in turn, directly influence the CO₂ dynamics across the MOF–polymer interface. This comprehensive analysis identifies critical aspects of interfacial pore structuring that control molecular diffusion, thereby paving the way for the rational design of MMMs with maximized performance.

Fig. 2 Workflow of our automated computational platform for constructing and analyzing MOF/polymer composites integrating DFT calculations for MOF slab model generation, FF-EMD for composites construction and gas translational and angular dynamics exploration, and FF-GCMC for gas adsorption study. Interfacial exploration includes structural characterization (atomic density profiles, radial distribution functions, Delaunay tessellation) and textural analysis (pore size distribution, void fraction, free pore mapping). Graph theory techniques (NetworkX) are employed to deliver comprehensive insights into the interfacial pore connectivity and possible molecular transport pathways (this new implementation is highlighted by the symbol*).

2. Automated computational platform for MMMs generation and analysis

In this section, we present an automated version of our computational suite, which integrates QM calculations at the Density Functional Theory (DFT) level, FF-EMD and FF-Monte Carlo (FF-MC) simulations, along with advanced post-simulation analyses incorporating graph theory techniques. This comprehensive platform is designed to: (i) construct atomistic models of MOF/polymer composites by combining MOF surface slabs with detailed atomistic polymer configurations; (ii) unravel the textural, structural, and morphological features, including pore connectivity, of the interfacial regions; and (iii) provide an in-depth characterization of gas adsorption and dynamics within these hybrid composites. Our workflow starts with DFT optimization of the MOF slab models cleaved from the periodic MOF structures. Herein, we applied this first step to generate MOF exposed surfaces perpendicularly to the direction of the MOF pore apertures, i.e. [001] (AlFFIVE-1-Ni), [100] CALF-20, [001] MIL-53-NH₂ and [110] Zr-fum-fcu. These MOF slab models are subsequently combined with an atomistic structure model of the 6FDA-DAM polymer generated via FF-EMD, with the objective of constructing MOF/polymer interfacial models, which are subsequently geometry-optimized using FF-EMD. Post-simulation structural and textural analysis of the MOFs/6FDA-DAM interfacial porosity includes atomic density profile, radial distribution functions (RDFs) between different MOF/polymer atom pairs, Delaunay tessellation, planarity assessments, pore size distribution (PSD), void fraction, and free pore mapping (see Fig. S1). Decisively, graph theory (using NetworkX)³⁴ is equally implemented to model interfacial pore topology, including interconnectivity, shape and dimension, providing an unprecedented topological analysis of the interfacial pore architecture. Adsorption properties can be further evaluated using FF-MC simulations, while FF-EMD calculations assess the CO₂ dynamics phenomena across the overall composite. Relative to our earlier workflow,¹³ the novelties here are: (i) a fully automated pipeline from atomistic polymer model generation, MOF/polymer composite generation to post-analysis, and (ii) a graph-theory module that enables an unprecedented quantification of the 3D interfacial pore-network topology (connectivity, dimensionality, and path lengths).

2.1. Construction of atomistic MOF surface models

The surface slab models for all the MOFs mentioned above are cleaved from their DFT-optimized crystal structures. Our focus is on MOF surfaces with pore entrances exposed to the polymer matrix. These exposed surface slab models of [001] ALFFIVE-1-Ni,²¹ [100] CALF-20 (ref. 28) and [110] Zr-fum-fcu²⁹ were adopted from our previous studies, ensuring consistency in methodology. These surface slab models are reminded in the ESI with their respective dimensions (see Fig. S2 to S4). The additional MIL-53-NH₂ surface slab model was cleaved along the [001] crystallographic plane of the DFT-geometry optimized MIL-53-NH₂ bulk structure (see Fig. S5). The terminations of these MOF surface models align with previously established models, ensuring a consistent framework for analysis.^17,21,28,29 Typically, the newly created [001] MIL-53-NH₂ slab model was terminated by –OH groups that mimics the scenario applied to Zr-fum-fcu.²⁹ All these MOF surface models were extended above 50 Å in length along the z direction, to prevent interactions between the surfaces in this direction. A vacuum box of z-length 20 Å was also added to avoid interactions between the MOF slab and its periodic images. We also ensured that the net dipole moment is zero in the z-direction for all these slab models by rebuilding the surface accordingly. Their geometries were further optimized using periodic DFT optimizations via the Quickstep module implemented in CP2K (Version 2023.1).^35,36 The Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional³⁷ was employed alongside Grimme's D3 dispersion correction,³⁸ which accounts for both short- and long-range interactions. Basis sets and pseudopotentials were sourced from the CP2K library, employing DZVP-MOLOPT-GTH for oxygen, carbon, aluminum, nickel, zinc, zirconium and fluorine atoms, and TZVP-MOLOPT-GTH for hydrogen atoms while GTH-PBE pseudopotentials were applied to all atoms. The grid settings featured a plane-wave cutoff of 500 Ry and a relative cutoff of 50 Ry. The SCF convergence criterion is set to 10⁻⁶, utilizing the Orbital Transformation (OT) method with DIIS minimization and a full kinetic preconditioner. Geometry optimizations were performed using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, with convergence criteria set to a maximum force of 2.0 × 10⁻⁴ hartree bohr⁻¹ and a maximum displacement of 2.0 × 10⁻⁴ bohr. The stress tensor was calculated analytically to ensure precise cell parameter optimization. These optimized slab models were systematically expanded to be in line with the dimensions of the AlFFIVE-1-Ni model. The final box sizes in the xy-plane were established as follows: 48.927 Å × 48.927 Å for AlFFIVE-1-Ni, 48.468 Å × 47.418 Å for CALF-20, 53.504 Å × 50.444 Å for Zr-fum-fcu, and 60.363 Å × 53.965 Å for MIL-53-NH₂.

The universal force-field (UFF) parameters were selected to describe the bonded and non-bonded terms of the MOF slab model. The use of the UFF was demonstrated previously to lead to reliable structure and adsorption properties for MOFs which do not contain any open metal as is the case for all MOFs investigated in this study as well as for previous MOF slab models we have explored including ALFFIVE-1-Ni²¹ and CALF-20.²⁸ The atomic partial charges for MIL-53-NH₂ surface model was computed using the Density Derived Electrostatic and Chemical (DDEC6) method, as implemented in the CHARGEMOL module.^39,40 following the same strategy we previously applied for the three previously reported MOF slab models, ALFFIVE-1-Ni,²¹ CALF-20,²⁸ and Zr-fum-fcu.²⁹ The full set of force field parameters and the charge values for each MOF are provided in the SI.

2.2. Construction of atomistic 6FDA-DAM polymer model

The 6FDA-DAM polymer model was developed using an automated in silico polymerization strategy that involves a construction of an all-atom monomer structure, initially built in Avogadro (version 1.2.0n)⁴¹ and subsequently geometry optimized using DFT calculations with the SCF/cc-pVDZ basis set. Atomic partial charges were calculated using the Restrained Electrostatic Potential (RESP) scheme, as implemented in Psi4 (version 1.4).⁴² The optimized monomer structure served as the basis for polymer chain construction. Chains consisting of 10 to 50 monomers were constructed, with terminations capped to mirror realistic chemical terminations (see Fig. S6). One chain end was capped by replacing a nitrogen atom with an oxygen atom, while the opposite end was terminated by attaching an –NH₂ group to a terminal carbon atom. The polymer chains were randomly packed into a cubic simulation box of 400 Å in size. A condensation protocol based on FF-EMD scheme was applied to transition the system into a liquid-like phase while preserving polydispersity and heterogeneity. The procedure involved an initial NVT ensemble FF-EMD equilibration at 300 K, followed by an NPT ensemble FF-EMD simulation at 600 K and 1 kbar to compress the system. Subsequently, the pressure gradually reduced to 1 bar under NPT/NVT ensembles FF-EMD conditions, completing the generation process (see Fig. S7). The polymer model was treated as fully flexible. Bonded interactions were modeled using the analytical expressions and parameters defined by the General Amber Force Field (GAFF),⁴³ while non-bonded interactions are described using a Lennard-Jones (LJ) potential combined with coulombic interactions. To perform structural relaxation and achieve realistic density profiles, the 6FDA-DAM polymer atomistic model underwent a series of 21FF-EMD equilibration steps.^13,23 This involved seven cycles of temperature and pressure adjustments: each cycle began with an NVT ensemble simulation at 600 K to promote chain mobility, followed by cooling to 300 K under NVT conditions to stabilize the polymer structure. Finally, an NPT ensemble simulation at 300 K was performed, with pressure ramping from 1 bar to 50 kbar in the initial cycles, then gradually decreasing to ambient pressure. Nosé–Hoover thermostats and barostats^44,45 with coupling constants of 0.1 ps were used for temperature and pressure regulation, and all simulations employed a 1 fs time step. The Lennard-Jones 12-6 cross-term were determined using the Lorentz–Berthelot mixing rule,⁴⁶ while electrostatic interactions were computed using the particle–particle particle-mesh solver (pppm),⁴⁷ with a tolerance level of 10⁻⁴. Both van der Waals and electrostatic interactions were modeled with a uniform cutoff distance of 12 Å. These simulations were conducted using the LAMMPS software package (version 2Aug2023).⁴⁸ The evolution of its density during the equilibration process were recorded to confirm proper structural relaxation (see Fig. S8). The resultant 6FDA-DAM atomistic model (1.17 g cm⁻³) was demonstrated to capture the experimental density (1.259 g cm⁻³)²³ and porosity of the corresponding polymer.

2.3. Construction of atomistic MOF/polymer composite models

The four distinct 6FDA-DAM/MOF composite models were constructed following the 21-cycles FF-EMD simulation equilibration protocol we described earlier^13,23 implying MD simulations performed in the NVT and NP_zT ensemble, where P_z is the pressure component in the direction normal to the MOF surface. The constructed composite models assemble about 40 wt% MOF loading in the 6FDA DAM polymer.²¹ To initiate the composite construction, the x- and y-dimensions of the 6FDA-DAM atomistic condensed model were adjusted to match the dimensions of the MOF surface model, and unwrapped along the z-direction to bring the polymer into contact with the MOF surface. During the NP_zT MD simulations, MOF acts as a piston on the polymer to compress and expand the system in z-direction. Both MOF and polymer components were treated as fully flexible, employing the force field parameters specified in the sections above. To ensure statistical reliability in the analysis, 3 independent replicas were generated for each composite system using the LAMMPS software package (version 2Aug2023).⁴⁸ Once the MOF/polymer composites were equilibrated analysis of the textural and structural features of their interfaces and the CO₂ translational and rotational were analyzed over the last 1 ns of the FF-EMD simulations.

2.4. Structural and textural analysis of the MOF/polymer interfacial porosity

An automated post-analysis suite was developed to assess key structural and textural descriptors of the interfacial pore network (see Fig. 2). This framework allows for effective determination of atomic density profiles, PSD, void fractions, and both 2D and 3D interfacial pore mappings. To define the 3D interfacial region, a novel approach was employed that uses Delaunay tessellation to delineate the boundaries between the MOF and the polymer phases. Surface atoms of each component are first identified, and a Voronoi diagram of the 3D structure is generated. The Delaunay tessellation discretizes the composite into three distinct regions: the MOF region, the polymer region, and the MOF/polymer interface. By tessellating the surfaces of both components, the method enables precise identification and structural characterization of the interfacial zone including the planarity deviation of the interfacial surfaces. Void fraction and pore network analyses were conducted using a uniform grid spacing of 0.6 Å in all three spatial dimensions. Grid points located beyond the sum of atomic radii plus an additional 0.8 Å were classified as free space, ensuring accurate identification of accessible void regions. Additionally, RDFs between MOF and polymer atom pairs were computed to provide detailed insights into interfacial interactions. All RDFs report normalize g(r) within the interfacial region defined by our Delaunay-based segmentation, where the MOF surface was discretized by selecting only the surface atoms of the framework.

2.5. Pore network topology analysis of MOF/polymer interfacial porosity

Graph theory techniques were implemented to gain deeper insights into the connectivity and topology of the interfacial pore networks defined from the Delaunay tessellation. Here the grid space was increased to 0.8 Å in all three directions in order to reduce the computational cost for the graph theory analysis without compromising the precise analysis of the interfacial pore structure. A probe particle with a radius of 1.3 Å was used to define accessible regions, effectively excluding spaces that are sterically inaccessible to small molecules. In the constructed graphs, nodes represent accessible probe particles, while edges are established between node pairs within 1.6 Å of each other. Edge weights were assigned based on the inverse of the distance (1/distance) to capture spatial proximity accurately. To analyze pore connectivity, community detection algorithms were included using the Clauset–Newman–Moore greedy modularity maximization method,⁴⁹ as implemented in NetworkX.³⁴ A resolution parameter of 0.1 is selected, favoring the formation of larger aggregates. Aggregates consisting of less than 10 nodes were excluded to focus on significant structural features. To characterize the interfacial pore networks, three representative graph theory network parameters were determined: (i) eccentricity,³⁴ redefined here as the network eccentricity normalized by the simulation box length along the x-axis, to capture the relative spatial extent of the pore network; (ii) assortativity^50,51 to evaluate the tendency of nodes to connect with similar or dissimilar nodes; and (iii) betweenness centrality^52–54 to identify plausible transport pathways for guest molecules.

2.6. Molecular gas adsorption and dynamics

FF-grand Canonical Monte Carlo (FF-GCMC) simulations were performed at 300 K and 10 bar gas pressure to load the four different MOF/6FDA-DAM composites with CO₂. These simulations were conducted using the Complex Adsorption and Diffusion Simulation Suite (CADSS) code.⁵⁵ In these calculations, the system undergoes 10⁷ equilibration steps followed by 10⁷ production steps to ensure accurate sampling of adsorption configurations. CO₂ molecules were modeled using the three-site linear rigid Lennard-Jones (LJ) TraPPE model.⁵⁶ Short-range interactions were truncated at a cutoff radius of 12 Å, while long-range electrostatic interactions were handled via the Ewald summation⁵⁷ technique to account for periodic boundary conditions accurately. The fugacity values required for GCMC simulations were derived from the Peng–Robinson equation of state.⁵⁸

CO₂ dynamics in these guest-loaded composites were further analyzed using FF-EMD simulations, which enabled the computation of velocity autocorrelation function (VACF) and the orientation distribution for the guest molecules providing insights into both translational and rotational dynamics of CO₂ in the MOF/polymer interfacial porosity. These interfacial analyses were performed with short FF-EMD production runs for 100 picoseconds (ps), with trajectory data collected every 10 fs for high-resolution temporal analysis. The CO₂ self-diffusivity (D^GK_s) was calculated using the Green–Kubo relation: , where v(t₀) and v(t) represent the reference and instantaneous velocities of CO₂ at times t₀ and t, respectively. The integrand, v(t₀)·v(t), is the VACF, which quantifies how the velocity of a particle at an initial time correlates with its velocity at a later time.^59–61 Rotational dynamics were assessed by analyzing the orientation distributions for each CO₂ molecule measuring the angle between its orientation vector at time t and at the reference time t₀. These orientation profiles were fitted using Gaussian functions, and the amplitude of angular deviation (θ_FWHM) was quantified throughout the determination of the Full Width at Half Maximum (FWHM) of the fit. The FWHM provides a direct measure of the angular dispersion, representing the extent to which the CO₂ molecules fluctuate about their average orientation. Larger FWHM values indicate greater rotational freedom, while narrower widths imply more constrained orientational alignment. By averaging the FWHM across all molecules in the interfacial region, the rotational landscape of CO₂ can be quantitatively compared between the different composites, offering insight into how the interfacial structuring impacts the rotational mobility of the guest molecules.

3. Results and discussion

3.1. Structural analysis of the MOF/6FDA-DAM composites

A structural analysis of the FF-EMD constructed atomistic models for the four different composites is provided in Fig. 3a–d. These representations include the density profiles for all composites plotted along the direction normal to the MOF surface, namely the z-axis, along with illustrations of the composites overlaid with Delaunay tessellation of the MOFs surfaces and plane representation of the polymer coverage.

Fig. 3 Structural analysis of MOF/6FDA-DAM composites. The atomic density profiles along the z-coordinates: (a) ALFFIVE-1-Ni/6FDA-DAM, (b) CALF-20/6FDA-DAM, (c) Zr-fum-fcu/6FDA-DAM, and (d) MIL-53-NH₂/6FDA-DAM composites. Insets show corresponding 2D and 3D pore mappings, highlighting the spatial distribution and connectivity of interfacial pores within each composite where the polymer surface is represented by the Green Delaunay tessellation. (e) The normalized radial distribution function, g(r), between the characteristic atoms of the MOF surfaces and the oxygen 6FDA-DAM atom pairs. Illustrative snapshots of the interfacial interactions for ALFFIVE-1-Ni/6FDA-DAM (f) and Zr-fum-fcu/6FDA-DAM (g), magenta lines represent the interactions between MOF surface atoms and the 6FDA-DAM carbonyl atoms.

The atomic density of 6FDA-DAM is shown to fluctuate around a mean value of 0.07 atoms per ų at both ends of the simulation box for all composites. In proximity to the MOF surfaces located at the center of the simulation box, the atomic density of 6FDA-DAM drops to zero. Interestingly, in the case of ALFFIVE-1-Ni showing the highest planarity (Fig. 1) the resulting composite does not exhibit any MOF/polymer overlap. 6FDA-DAM distributes homogeneously at the ALFFIVE-1-Ni surface, this geometry being held by means of continued moderate van der Waals interactions between its terminal F-atom and the carbonyl functional group (C=O) of the 6FDA-DAM polymer with separating distance over 3.5 Å as shown in the corresponding RDF plotted in Fig. 3e. Fig. 3f delivers an illustration of the closest contacts between the fluorine atoms of ALFFIVE-1-Ni and the carbonyl groups of 6FDA-DAM. On the other hand, for the other 3MOF slab models showing a substantial planarity deviation, the two components overlap at the interface accompanied with interactions between the –OH terminal groups of the MOF and the carbonyl groups of 6FDA-DAM as evidenced in the RDF plots (Fig. 3e) and their representative snapshots (Fig. 3f). The MOF/polymer overlap length in the 2D atomic density profile is slightly longer for the MOF surface with the higher planarity deviation: Zr-fum-fcu/6FDA-DAM (8.6 Å) > MIL-53-NH₂/6FDA-DAM (7.5 Å) > CALF-20/6FDA-DAM (3.3 Å).

When analyzing the surface tessellation, we observed that the MIL-53-NH₂ and Zr-fum-fcu surfaces exhibit a tessellation pattern similar to that of 6FDA-DAM, indicating that the polymer effectively adapts to the shape of the MOF surface which deviates with the behavior observed for CALF-20 and even more for ALFFIVE-1-Ni for which the surface tessellations are more erratic. To further characterize these different behaviors, we computed the 3D interface gap by determining the Delaunay tessellation surface area and dividing it by the free space volume between the 6FDA-DAM and the MOF Delaunay surfaces. The calculated interface gaps for each composite are 4.1 ± 0.05 Å (CALF-20), 3.8 ± 0.4 Å (ALFFIVE-1-Ni), 3.0 ± 0.1 Å (Zr-fum-fcu), and 2.8 ± 0.1 Å (MIL-53-NH₂). This newly derived descriptor serves as a reliable metrics to assess the overall 3D polymer adhesion to the MOF surface, with MIL-53-NH₂ and Zr-fum-fcu exhibiting the shortest interfacial gaps in line with a better adaptability of the polymer conformation to the morphology of their surfaces.

3.2. Textural analysis of the MOF/6FDA-DAM interfaces

The interfacial PSD plots depicted in Fig. 4a for the different MOF/6FDA-DAM composites evidence that the CALF-20/6FDA-DAM composite exhibits the largest pore size contribution with a dominant peak around 6 Å. This suggests that the resulting interface maintains relatively larger voids compared to the scenario obtained for the other MOFs. The ALFFIVE-1-Ni/6FDA-DAM interface presents a distinct peak at approximately 4.3 Å, with a sharper decline toward smaller pore sizes, indicating a more homogeneous distribution of the interfacial pore sizes compared to CALF-20. The PSDs of the MIL-53-NH₂/6FDA-DAM and Zr-fum-fcu/6FDA-DAM interfaces show a broader distribution with significant contributions across a range of pore sizes, peaking between 3 and 5 Å. This indicates a more heterogeneous interfacial pore size distribution compared to the other composites.

Fig. 4 Textural analysis of the CALF-20/6FDA-DAM (green), ALFFIVE-1-Ni/6FDA-DAM (blue), Zr-fum-fcu/6FDA-DAM (orange) and, MIL-53-NH₂/6FDA-DAM (cyan) composites. (a) Interfacial PSD plots. (b) Void fraction profiles as a function of the distance from the MOF surface, plotted along the z-direction of the composites. The position z = 0 corresponds the metal atom closest to the MOF surface. (c) Representative illustration of the free porosity over the xy plane for different z values when one moves away from the MOF surface, where colored regions represent different void domains, and white regions denote the presence of polymers.

Fig. 4b illustrates the void fraction distribution of the composites in the (x,y) plane for different z-values considered along the direction normal to the MOF surface, for the different MOF/6FDA-DAM interfaces, providing further confirmation of the interfacial porosity trends observed in the PSD analysis. These plots reveal how the interfacial free volume is distributed when we move away from the metal atoms closest to the MOF surface; Ni (ALFFIVE-1-Ni), Zn (CALF-20), Zr (Zr-fum-fcu) and Al (MIL-53-NH₂) respectively. The CALF-20/6FDA-DAM interface exhibits the largest interfacial void fraction, maintaining values above 0.4 and peaking beyond 0.6. This trend consistent with its associated broad PSD profile, highlights that CALF-20 surface favours the formation of an interface with substantial free pore volume. The ALFFIVE-1-Ni/6FDA-DAM interface displays a pronounced peak in void fraction around 4 Å, aligning with its narrow PSD peak at the same distance. This suggests that the interface is well-structured with a defined porosity region as we suggested in earlier works.^21,22 The MIL-53-NH₂/6FDA-DAM and Zr-fum-fcu/6FDA-DAM interfaces exhibit a lower and more evenly distributed void fraction, which aligns with their broader porosity depicted in the PSD profile. This indicates that these two interfaces are somehow less porous with a rather heterogeneous distribution of interfacial small-to-large pores.

To visualize the spatial distribution of interfacial pores, Fig. 4c provides representative (x,y) plane free pore mappings drawn at different z values for the different MOF/6FDA-DAM composites. Additional maps from independent simulation replicas are provided in Fig. S9–S12. These mappings illustrate how free volume is spatially distributed at the interface, further complementing the PSD and void fraction analyses. Each row corresponds to a distinct MOF/6FDA-DAM composite, with different colours representing individual pore domains, this distinction being discretized using the function Clauset–Newman–Moore⁴⁹ greedy modularity maximization with a resolution of 0.1. The CALF-20/6FDA-DAM interface exhibits rather large and well-connected pore domains, in line with its high void fraction and broad PSD distribution. At the opposite, MIL-53-NH₂/6FDA-DAM interface shows the lowest pore interconnectivity with a fragmented distribution of small to medium-sized pores. Intermediate profiles are observed for the two other ALFFIVE-1-Ni/6FDA-DAM and Zr-fum-fcu/6FDA-DAM interfaces characterized by the co-existence of pore domains with rather large connectivity and isolated small pores.

3.3. Interfacial pore network topology analysis of the MOF/6FDA-DAM interfaces

Graph theory analysis enabled to characterize in-depth the interfacial pore network interconnectivity for all MOF/6FDA-DAM composites with the evaluation of a series of representative network parameters, e.g. assortativity (how asymmetric the pore domain is), betweenness centrality (how concentrated the pore domain is) and eccentricity (how large the pore domain is) (Fig. 5). Each porous domain identified in the MOF/6FDA-DAM interfaces was analysed independently and weighted as a function of its volume to obtain the closest representation of the pore network. Typically, the CALF-20/6FDA-DAM interface which contains large interconnected interfacial pore domains with small pores being mostly absent results in the highest calculated assortativity and eccentricity values. At the opposite the fragmented pore network of MIL-53-NH₂/6FDA-DAM leads to the lowest degree of connectivity, reflected by its resulting high betweenness centrality and low eccentricity values. ALFFIVE-1-Ni/6FDA-DAM interface shows less connected but still highly structured pore domains with moderate betweenness centrality, assortativity and eccentricity values. Similar metrics are obtained for the Zr-fum-fcu/6FDA-DAM interface in line with equivalent free pore mappings reported in Fig. 4 for these two interfaces. These detailed analyses emphasize the critical role of MOF surface roughness in shaping interfacial pore network connectivity.

Fig. 5 3D visualizations of the interfacial pore networks for CALF-20 (green), ALFFIVE-1-Ni (blue), MIL-53-NH₂ (cyan) and Zr-fum-fcu (yellow) composites, highlighting differences in pore connectivity and distribution. Radar plots summarize three key graph theory network parameters: assortativity, eccentricity, and betweenness centrality.

3.4. Interfacial CO₂ dynamics

We next examined the influence of pore network topology on interfacial CO₂ dynamics. Quantitative metrics derived from the simulations, i.e. the translational mobility from the velocity autocorrelation functions (VACFs) and the rotational freedom from the angular distributions are reported in Fig. S13 and S14, respectively. Possible correlations between the representative interfacial pore network parameters and the CO₂ dynamics were evaluated and summarized in Fig. 6. Our analysis revealed that eccentricity (and assortativity in Fig. S15) exhibits a positive linear correlation with the CO₂ translational self-diffusivity D^GK_s (Fig. 6a) whereas betweenness centrality displays an inverse relationship (Fig. 6b). The CALF-20/6FDA-DAM interface enables an effective CO₂ interfacial motion in line with its large interfacial interconnected pore domains characterized by high assortativity and eccentricity. Moderate CO₂ interfacial dynamics were observed in both the ALFFIVE-1-Ni/6FDA-DAM and Zr-fum-fcu/6FDA-DA composites, reflecting the limited interconnectivity of their pore domains. Their corresponding moderate betweenness centrality and eccentricity values indicate that, although migration pathways are present, they are more constrained than in CALF-20/6FDA-DAM leading to a more limited CO₂ translational dynamics. Finally, the MIL-53-NH₂/6FDA-DAM interface exhibits the slowest CO₂ translational mobility, in line with its fragmented and highly dispersed interfacial pore network. This suggests that CO₂ molecules face larger mobility constraints, further exacerbated by the presence of isolated, small pores that hinder guest dynamics. Notably, the rotational freedom of interfacial CO₂ molecules correlates positively with D^GK_s, indicating that reduced interfacial steric hindrance enhances both rotational and translational dynamics at the interface (Fig. 6c). This overall analysis reveals that more interconnected the interfacial pore network, faster is the CO₂ dynamics (case of CALF-20/6FDA-DAM) while highly fragmented porosity, leads to most restricted interfacial CO₂ translational and rotational mobility (case of MIL-53-NH₂/6FDA-DAM).

Fig. 6 Correlation between the interfacial CO₂ self-diffusivity (from the Green–Kubo relation) (D^GK_s) and (a) eccentricity, (b) betweenness centrality, and (c) CO₂ rotational dynamics quantified by θ_FWHM. Error bars on the graph theory network parameters represent their standard deviations evaluated over the 3 interfacial models constructed for all composites. The CO₂ dynamics metrics are equally averaged over the 3 interfacial models. Dashed lines are included to guide the eye.

Our previous work demonstrated that an accelerated CO₂ interfacial migration in MOF-polymer MMM likely boosts the molecular transport all along the MMM and hence a faster CO₂ permeability.²¹ Accounting for differences in MOF loading in the MMM, working gas pressure and temperature conditions employed, previously reported experimental CO₂ permeability data positions (001) ALFFIVE-1-Ni/6FDA-DAM achieves ∼1295.7 barrer at 60 wt% (2 bar, 35 °C),²² above MIL-53-NH₂/6FDA-DAM (∼650 barrer at 20 wt%, 3 bar, 25 °C),⁶² This aligns with our 3D interfacial-network analysis, where higher connectivity (higher eccentricity, lower betweenness) yields faster interfacial CO₂ dynamics for ALFFIVE-1-Ni/6FDA-DAM and suggest based on our discussion above a qualitative trend for CO₂ permeability as follows: MIL-53-NH₂/6FDA-DAM < Zr-fum-fcu/6FDA-DAM ≲ ALFFIVE-1-Ni/6FDA-DAM < CALF-20/6FDA-DAM.

4. Conclusions

This study provides a comprehensive molecular understanding of how MOF surface morphology governs the topology of interfacial pore networks in MOF/polymer MMMs. By developing an automated, graph theory-enhanced molecular simulation platform capable of high-resolution analysis of complex interfacial architectures, we demonstrate that MOF surface planarity and roughness play a decisive role in shaping the dimensionality, connectivity, and spatial organization of interfacial pores. Importantly, we show that highly interconnected and continuous interfacial pore networks enable more efficient interfacial translational and rotational CO₂ dynamics, whereas fragmented structures severely impede guest mobility. Beyond delivering unprecedented microscopic insights, this work introduces a new paradigm: the deliberate tuning of MOF surface morphology, typically via surface grafting, as a powerful strategy to engineer interfacial pore networks and optimize interfacial gas dynamics. We note that interfacial CO₂ dynamics cannot be measured directly; experimental validation presently relies on overall permeability trends. While a rigorous, system-by-system validation is outside our scope here, our mechanistic link that more-interconnected interfacial networks promote higher interfacial mobility and expected faster CO₂ permeability in the overall MMM.²¹ By unlocking this largely unexplored strategy, our findings open new avenues for the rational design of next-generation MMMs, offering opportunities for advancing high-performance and energy-efficient separation technologies.

Author contributions

A. D. M. and G. M. conceived the idea of the project, designed the simulation strategy and wrote the first draft of the manuscript. D. F. constructed and analyzed the MOF/polymer composite atomistic models. S. N. performed the polymer construction simulations. G. M. supervised the computational work and M. E. participated in the discussion. All authors contributed to the revision of the manuscript.

Conflicts of interest

The authors have no conflict of interest to declare.

Data availability

Code availability: The central codes used in this paper are VASP,⁶³ CP2K,³⁵ LAMMPS,⁴⁸ CADSS,⁵⁵ and polymatic.⁶⁴ Detailed information related to the license and user guide are available in the referenced papers and on their websites. We have deposited all input files used in this work (LAMMPS scripts, UFF parameter files, and partial charges), together with post-processing scripts (graph-theory analysis), in the repository https://github.com/AlejandroDiazMarquez/AutoMMM. The automated polymerization scripts are in the repository https://github.com/supriyonaskar/Polymer.git.

All data needed to evaluate the conclusions in the paper are present in the paper and/or the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sc04241k.

Acknowledgements

The research reported in this publication was supported by King Abdullah University of Science and Technology (KAUST). The computational work reported in this publication was supported by KAUST (CCF 1972 project) and performed using HPC resources from GENCI-CINES (Grant A0180907613). GM thanks Institut Universitaire de France for the Senior Chair.

References

1. S. D. Sholl and L. P. Ryan, Seven chemical separations to change the world, Nature, 2016, 532, 435–437.
2. S. J. Miller, W. J. Koros and D. Q. Vu Mixed matrix membrane technology: enhancing gas separations with polymer/molecular sieve composites. in Studies in Surface Science and Catalysis, ed. Xu, R., Gao, Z., Chen, J. and Yan, W., Elsevier, 2007, vol. 170, pp. 1590–1596.
3. B. D. Freeman, Basis of permeability/selectivity tradeoff relations in polymeric gas separation membranes, Macromolecules, 1999, 32, 375–380.
4. L. M. Robeson, The upper bound revisited, J. Membr. Sci., 2008, 320, 390–400.
5. H. B. Park, J. Kamcev, L. M. Robeson, M. Elimelech and B. D. Freeman, Maximizing the right stuff: The trade-off between membrane permeability and selectivity, Science, 2017, 356, 1138–1148.
6. R. Lin, B. Villacorta Hernandez, L. Ge and Z. Zhu, Metal organic framework based mixed matrix membranes: An overview on filler/polymer interfaces, J. Mater. Chem. A, 2018, 6, 293–312.
7. X. Guo, Z. Qiao, D. Liu and C. Zhong, Mixed-matrix membranes for CO2 separation: Role of the third component, J. Mater. Chem. A, 2019, 7, 24738–24759.
8. N. Tien-Binh, H. Vinh-Thang, X. Y. Chen, D. Rodrigue and S. Kaliaguine, Polymer functionalization to enhance interface quality of mixed matrix membranes for high CO2/CH4 gas separation, J. Mater. Chem. A, 2015, 3, 15202–15213.
9. G. Yu, et al., Constructing Connected Paths between UiO-66 and PIM-1 to Improve Membrane CO 2 Separation with Crystal-Like Gas Selectivity, Adv. Mater., 2019, 31, 1806853.
10. H. Molavi and A. Shojaei, Mixed-Matrix Composite Membranes Based on UiO-66-Derived MOFs for CO 2 Separation, ACS Appl. Mater. Interfaces, 2019, 11, 9448–9461.
11. C. Ma and J. J. Urban, Hydrogen-Bonded Polyimide/Metal-Organic Framework Hybrid Membranes for Ultrafast Separations of Multiple Gas Pairs, Adv. Funct. Mater., 2019, 29, 1903243.
12. Y. Cheng, et al., Advances in metal–organic framework-based membranes, Chem. Soc. Rev., 2022, 51, 8300–8350.
13. R. Semino, N. A. Ramsahye, A. Ghoufi and G. Maurin, Microscopic Model of the Metal-Organic Framework/Polymer Interface: A First Step toward Understanding the Compatibility in Mixed Matrix Membranes, ACS Appl. Mater. Interfaces, 2016, 8, 809–819.
14. Q. Qian, et al., Mixed-Matrix Membranes Formed from Imide-Functionalized UiO-66-NH2 for Improved Interfacial Compatibility, ACS Appl. Mater. Interfaces, 2019, 11, 31257–31269.
15. R. Semino, N. A. Ramsahye, A. Ghoufi and G. Maurin, Role of MOF surface defects on the microscopic structure of MOF/polymer interfaces: A computational study of the ZIF-8/PIMs systems, Microporous Mesoporous Mater., 2017, 254, 184–191.
16. D. Fan, et al., Engineering MOF surface defects in mixed matrix membranes: An effective strategy to enhance MOF/polymer adhesion and control interfacial gas transport, J. Membr. Sci. Lett., 2022, 2, 100029.
17. R. Semino, J. C. Moreton, N. A. Ramsahye, S. M. Cohen and G. Maurin, Understanding the origins of metal-organic framework/polymer compatibility, Chem. Sci., 2018, 9, 315–324.
18. D. Fan, et al., Is Porosity at the MOF/Polymer Interface Necessarily an Obstacle to Optimal Gas-Separation Performances in Mixed Matrix Membranes?, ACS Mater. Lett., 2021, 3, 344–350.
19. S. R. Tavares, et al., Computationally Assisted Assessment of the Metal-Organic Framework/Polymer Compatibility in Composites Integrating a Rigid Polymer, Adv. Theory Simul., 2019, 2, 1900116.
20. I. D. Carja, et al., Insights into the Enhancement of MOF/Polymer Adhesion in Mixed-Matrix Membranes via Polymer Functionalization, ACS Appl. Mater. Interfaces, 2021, 13, 29041–29047.
21. A. Ozcan, et al., Tuning MOF/polymer interfacial pore geometry in mixed matrix membrane for upgrading CO 2 separation performance, Sci. Adv., 2024, 10, 5846.
22. S. Jit Datta, et al., Rational design of mixed-matrix metal-organic framework membranes for molecular separations, Science, 2022, 376, 1080–1087.
23. M. Z. Ahmad, et al., Enhanced gas separation performance of 6FDA-DAM based mixed matrix membranes by incorporating MOF UiO-66 and its derivatives, J. Membr. Sci., 2018, 558, 64–77.
24. E. B. Doǧan, G. Maurin and M. G. Ahunbay, Atomistic Insight into the Interfacial Structuring of ZIF-67 MOF/Polymer Composites and Their Propylene–Propane Adsorption Properties, J. Phys. Chem. C, 2023, 127, 20491–20502.
25. R. Hou, et al., Highly permeable and selective mixed-matrix membranes for hydrogen separation containing PAF-1, J. Mater. Chem. A, 2020, 8, 14713–14720.
26. J. Zhao, et al., Manipulating the interfacial interactions of composite membranes via a mussel-inspired approach for enhanced separation selectivity, J. Mater. Chem. A, 2015, 3, 19980–19988.
27. D. Fan, S. Naskar and G. Maurin, Unconventional mechanical and thermal behaviours of MOF CALF-20, Nat. Commun., 2024, 15, 3251.
28. S. Naskar, D. Fan, A. Ghoufi and G. Maurin, Microscopic insight into the shaping of MOFs and its impact on CO2 capture performance, Chem. Sci., 2023, 14, 10435–10445.
29. R. Hardian, et al., Design of Mixed-Matrix MOF Membranes with Asymmetric Filler Density and Intrinsic MOF/Polymer Compatibility for Enhanced Molecular Sieving, Adv. Mater., 2024, 36, 2314206.
30. T. H. Lee, et al., Disclosing the Role of Defect-Engineered Metal–Organic Frameworks in Mixed Matrix Membranes for Efficient CO2 Separation: A Joint Experimental-Computational Exploration, Adv. Funct. Mater., 2021, 31, 2103973.
31. P. Serra-Crespo, et al., Interplay of metal node and amine functionality in NH2-MIL-53: Modulating breathing behavior through intra-framework interactions, Langmuir, 2012, 28, 12916–12922.
32. A. Pustovarenko, et al., Nanosheets of Nonlayered Aluminum Metal–Organic Frameworks through a Surfactant-Assisted Method, Adv. Mater., 2018, 30, 1707234.
33. C. B. Barber, D. P. Dobkin and H. Huhdanpaa, The quickhull algorithm for convex hulls, ACM Trans. Math. Software, 1996, 22, 469–483.
34. A. A. Hagberg, D. A. Schult and P. J. Swart, Exploring Network Structure, Dynamics, and Function using NetworkX, in Proc. 7th Python in Science Conf., ed. G. Varoquaux, T. Vaught and J. Millman, Pasadena, CA USA, 2008, pp. 11–15.
35. T. D. Kühne, et al., CP2K: An electronic structure and molecular dynamics software package -Quickstep: Efficient and accurate electronic structure calculations, J. Chem. Phys., 2020, 152, 070401.
36. J. Vandevondele, et al., Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach, Comput. Phys. Commun., 2005, 167, 103–128.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem., 2006, 27, 1787–1799.
39. T. A. Manz and D. S. Sholl, Improved atoms-in-molecule charge partitioning functional for simultaneously reproducing the electrostatic potential and chemical states in periodic and nonperiodic materials, J. Chem. Theory Comput., 2012, 8, 2844–2867.
40. T. A. Manz and N. G. Limas, Introducing DDEC6 atomic population analysis: Part 1. Charge partitioning theory and methodology, RSC Adv., 2016, 6, 47771–47801.
41. M. D. Hanwell, et al., Avogadro: an advanced semantic chemical editor, visualization, and analysis platform, J. Cheminf., 2012, 4, 17.
42. D. G. A. Smith, et al., PSI4 1.4: Open-source software for high-throughput quantum chemistry, J. Chem. Phys., 2020, 152, 184108.
43. K. G. Sprenger, V. W. Jaeger and J. Pfaendtner, The general AMBER force field (GAFF) can accurately predict thermodynamic and transport properties of many ionic liquids, J. Phys. Chem. B, 2015, 119, 5882–5895.
44. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519.
45. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A, 1985, 31(3), 1695–1697.
46. H. A. Lorentz, Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase, Ann. Phys., 1881, 248, 127–136.
47. R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, CRC Press, Boca Raton, 2021.
48. A. P. Thompson, et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
49. A. Clauset, M. E. J. Newman and C. Moore, Finding community structure in very large networks, Phys. Rev. E, 2004, 70, 66111.
50. M. E. J. Newman, Mixing patterns in networks, Phys. Rev. E, 2003, 67, 26126.
51. J. G. Foster, D. V. Foster, P. Grassberger and M. Paczuski, Edge direction and the structure of networks, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 10815–10820.
52. U. Brandes and C. Pich, Centrality estimation in large networks, Int. J. Bifurcation Chaos, 2007, 17, 2303–2318.
53. U. Brandes, On variants of shortest-path betweenness centrality and their generic computation, Soc. Netw., 2008, 30, 136–145.
54. U. Brandes, A faster algorithm for betweenness centrality, J. Math. Sociol., 2001, 25, 163–177.
55. Q. Yang and C. Zhong, Molecular Simulation of Carbon Dioxide/Methane/Hydrogen Mixture Adsorption in Metal–Organic Frameworks, J. Phys. Chem. B, 2006, 110, 17776–17783.
56. J. J. Potoff and J. I. Siepmann, Vapor-liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen, AIChE J., 2001, 47, 1676–1682.
57. P. P. Ewald, Die Berechnung optischer und elektrostatischer Gitterpotentiale, Ann. Phys., 1921, 369, 253–287.
58. B. Zohuri Properties of Pure Substances, in Physics of Cryogenics, ed. B. Zohuri, Elsevier, 2018, pp. 53–79,  DOI:10.1016/B978-0-12-814519-7.00002-1.
59. M. S. Green, Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids, J. Chem. Phys., 1954, 22, 398–413.
60. R. Kubo, Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems, J. Phys. Soc. Jpn., 1957, 12, 570–586.
61. B. J. Alder and T. E. Wainwright, Decay of the Velocity Autocorrelation Function, Phys. Rev. A, 1970, 1, 18–21.
62. A. Sabetghadam, et al., Metal Organic Framework Crystals in Mixed-Matrix Membranes: Impact of the Filler Morphology on the Gas Separation Performance, Adv. Funct. Mater., 2016, 26, 3154–3163.
63. J. Hafner, Ab-initio simulations of materials using VASP: Density-functional theory and beyond, J. Comput. Chem., 2008, 29, 2044–2078.
64. L. J. Abbott, K. E. Hart and C. M. Colina, Polymatic: A generalized simulated polymerization algorithm for amorphous polymers, Theor. Chem. Acc., 2013, 132, 1–19.

---