Modeling of Gas Transport through Polymer/MOF Interfaces: A Microsecond-Scale Concentration Gradient-Driven Molecular Dynamics Study.

Membrane-based separation technologies offer a cost-effective alternative to many energy-intensive gas separation processes, such as distillation. Mixed matrix membranes (MMMs) composed of polymers and metal-organic frameworks (MOFs) have attracted a great deal of attention for being promising systems to manufacture durable and highly selective membranes with high gas fluxes and high selectivities. Therefore, understanding gas transport through these MMMs is of significant importance. There has been longstanding speculation that the gas diffusion behavior at the interface formed between the polymer matrix and MOF particles would strongly affect the global performance of the MMMs due to the potential presence of nonselective voids or other defects. To shed more light on this paradigm, we have performed microsecond long concentration gradient-driven molecular dynamics (CGD-MD) simulations that deliver an unprecedented microscopic picture of the transport of H2 and CH4 as single components and as a mixture in all regions of the PIM-1/ZIF-8 membrane, including the polymer/MOF interface. The fluxes of the permeating gases are computed and the impact of the polymer/MOF interface on the H2/CH4 permselectivity of the composite membrane is clearly revealed. Specifically, we show that the poor compatibility between PIM-1 and ZIF-8, which manifests itself by the presence of nonselective void spaces at their interface, results in a decrease of the H2/CH4 permselectivity for the corresponding composite membrane as compared to the performances simulated for PIM-1 and ZIF-8 individually. We demonstrate that CGD-MD simulations based on an accurate atomistic description of the polymer/MOF composite is a powerful tool for characterization and understanding of gas transport and separation mechanisms in MMMs.

Introduction

Membrane technology plays an important role in today’s industrial gas separation processes and has paramount economic importance.[1−4] The efficiency of membrane-based separation technologies reduces the cost and environmental footprint of many industrial processes.[5] The annual size of the membrane-based gas separations market, which was in the range of US$1–1.5 billion in 2017, is a concrete example of their impact.[6] Some examples of major membrane-based gas separation processes include hydrogen recovery from various off-gas streams, on-site nitrogen separation from air, natural gas sweetening, and olefin recovery from nitrogen-containing petrochemical vent gas streams.[7] In addition to the gas separation market, membrane-based applications in water desalination, organic solvent nanofiltration and waste water treatment are also attracting a great deal of attention.[8]

There has been substantial progress in the development of membranes for gas separations over the past 2 decades;[9−16] however, there are still a number of long-standing problems. The main challenge is to overcome the trade-off between permeability and selectivity, which has been illustrated by Robeson’s upper bound.[17,18] Even though polymeric membranes dominate the majority of current membrane-based applications, they are bound by this permeability–selectivity trade-off. Combining polymers and metal–organic frameworks (MOFs) in the form of a mixed matrix membrane (MMM) has been proposed as an alternative strategy, and this has shown significant promise.[19−23] This approach aims to take advantage of both the good processability of polymers and the excellent separation performances of crystalline porous MOFs. On the other hand, gas transport dynamics at the polymer/MOF interface is expected to play a determining role in the performance of the composite membrane and understanding the effects of polymer/MOF compatibility on gas separation is not a trivial task.[24] Molecular simulations can be used to quantify and characterize such interface effects in polymer/MOF composites provided that (i) accurate methods for the computation of the flux of permeants are employed and (ii) realistic atomistic models of the polymer/MOF interface are constructed.

The earliest molecular simulation study of gas separation in a polymer/MOF MMM was reported by Zhang et al.[25] on H2/CO2 in polybenzimidazole (PBI)/ZIF-7 MMMs using equilibrium molecular dynamics simulations. Velioglu et al.[26] and Altintas et al.[27] recently reported the separations of H2/CH4 and CO2/CH4 mixtures in polymer/MOF MMMs by carrying out screening calculations based on molecular simulations. However, these computational studies predicted the separation performance of the polymer/MOF composites based on the individual constituents of the MMMs by assuming ideal polymer/MOF compatibility and did not consider the impact of the interface on the transport properties. Furthermore, there are macroscopic models of permeation widely used to predict the permeability of MMMs based on the permeability data available for their constituent materials. Their applicability and limitations are discussed elaborately in a comprehensive review by Vinh-Thang et al.[28] These models are usually based on analogies with continuum models that define the thermal or electrical conduction in a heterogeneous medium, such as the so-called serial resistance model, but they usually do not take into account the interface effects.

Accurate atomistic polymer/MOF interface models have been developed by Semino et al.[29] and first applied to the PIM-1/ZIF-8 MMM to investigate the surface compatibility (i.e., the affinity) between the polymer and the MOF. These models have also been successfully applied to other polymer/MOF pairs, such as poly(vinyl alcohol)/HKUST-1[30] and 6-FDA-DAM/UiO-66 MMMs.[31] Further, these models were considered in an investigation of the transient concentration of CO2 in 6-FDA-DAM/ZIF-8 MMMs[32] by molecular modeling and IR microimaging. Despite these efforts to understand the molecular basis of polymer/MOF compatibility, no correlation has yet been found between compatibility and the performance of these composites for different applications. Even though the presence of microscopic sized voids at the interface has been clearly related to the formation of brittle membranes (i.e., nonoptimal mechanical properties),[29] the speculation that the presence of these voids might lead to a reduction in selectivity has never been confirmed. Conversely, the absence of these voids is a signature of good polymer/MOF compatibility, but this is not necessarily related to the good performance of the corresponding membrane for a particular application.

To address this still open question, here we report concentration gradient-driven molecular dynamics (CGD-MD) simulations of H2 and CH4 transport through a realistic atomistic model of the polymer/MOF membrane with a specific focus on gas transport properties through the interfaces as well as along the individual components of the MMM. CGD-MD is a nonequilibrium MD method recently developed to study the transport and separation of fluids through membranes.[33] The advantages of employing the CGD-MD method for the separation of gas mixtures over equilibrium MD approaches have been recently demonstrated.[34] However, it has long been speculated that the decrease in the separation of performances of MMMs might be related to the presence of defects at the interface between the two components. Our work provides a first clear confirmation at the microscopic level that this is really the case and we can equally quantify the negative impact of a poor compatibility at the interface.

As a proof of concept, the PIM-1/ZIF-8 composite was considered as a model membrane for H2/CH4 separation. PIM-1 (polymer of intrinsic microporosity-1) is a member of a group of microporous glassy polymers introduced by McKeown et al.[35] They are rigid, highly contorted spirobisindane-based ladder polymers, and their backbones have essentially no rotational freedom. This results in relatively large Brunauer–Emmett–Teller (BET) areas (∼800 m2/g)[36] and high permanent gas permeabilities. ZIF-8 is one the most studied MOF material and is known to have exceptional thermal and chemical stabilities.[36−38] It has large cages of 11.6 Å connected by 3.4 Å pore apertures and has been applied for various gas separation processes.[39−43] Furthermore, H2/CH4 separation by membranes is part of the $200 million/year hydrogen recovery market, which is substantially dominated by polysulfone and polyimide membranes.[6] Due to the relative difference in the size of H2 and CH4 molecules (kinetic diameters of 2.8 and 3.8 Å, respectively), these molecules are expected to exhibit distinct transport properties in the different regions of the polymer/MOF MMM, including the interfaces. We demonstrate below that this is indeed the case in the PIM-1/ZIF-8 membrane.

Models and Computational Details

Construction of the ZIF-8, PIM-1, and Composite PIM-1/ZIF-8 Membrane Models

The ZIF-8 membrane was derived from a previous work.[29] It consists of a [011] surface, terminated by −OH and −H groups, as per the dissociative adsorption of water, the standard solvent considered in the ZIF-8 synthesis, on the under-coordinated sites. This model was optimized at the density functional theory level, and it is periodic in the x- and y-directions. The net dipole in the z-direction (i.e., the direction normal to the membrane) is zero. The dimensions of the ZIF-8 model are 5.0, 4.8, and 9.8 nm in the x-, y-, and z-directions, respectively.

Different approaches were reported in the literature for the generation of polymer models.[44−46] Here, the construction of the PIM-1 membrane was performed using the in silico polymerization approach developed by Abbott et al. as implemented in the polymatic code,[47] which was previously employed to build different polymer models, including PIM-1.[29,47−49] The length of the resulting PIM-1 membrane was 20.8 nm in the z-direction. The composite PIM-1/ZIF-8 membrane was further constructed by putting together the models of ZIF-8 and PIM-1 in a simulation box and letting the polymer equilibrate in the presence of the MOF. This was achieved by a series of MD simulations, including annealing steps and a rapid compression followed by a slow decompression. Further details of this procedure can be found elsewhere.[29] The polymer/MOF model was further unwrapped in the z-direction and the polymer slab was duplicated on each side of the MOF in such a way that the MOF was located between two polymer slabs. The resulting MMM of 52.4 nm in the z-direction was further equilibrated by MD simulations. The constructed three membranes (ZIF-8, PIM-1, and composite PIM-1/ZIF-8) are illustrated in Figure .

Figure 1

Illustration of the (a) ZIF-8 membrane, (b) PIM-1 membrane, and (c) composite PIM-1/ZIF-8 membrane structural models used in the CGD–MD simulations. Color code: C (gray), O (red), N (blue), Zn (steel blue), and H (white).

Illustration of the (a) ZIF-8 membrane, (b) n class="Gene">PIM-1 membrane, and (c) composite PIM-1/ZIF-8 membrane structural models used in the CGD–MD simulations. Color code: C (gray), O (red), N (blue), Zn (steel blue), and H (white).

Modeling of the Gas Transport

Simulations of gas transport through the membranes were performed using GROMACS-5.1.2 simulation package[50] patched with a modified version of the PLUMED-2 enhanced sampling plug-in[51] to enable running the CGD-MD simulations.[33] In CGD-MD simulations, a concentration gradient between the feed and the permeate sides is created, which facilitates the transport of molecules across the membrane. The molecular fluxes can then be directly calculated from the CGD-MD simulations. To generate the concentration gradient across the membrane, the density of fluid molecules within designated volumes located at the inlet and outlet of the membrane are taken as collective variables and maintained at a target value with an external biasing scheme. An illustration of the CGD-MD setup and the parameters are provided in Figure S1 and Table S1, respectively. Further details of the method can be found elsewhere.[33,52] In all CGD-MD simulations, the membranes were placed in the middle of the simulation box and void space was added to both sides of the membranes, resulting in simulation box lengths of 40.8 nm for the PIM-1 membrane, 29.8 nm for the ZIF-8 membrane, and 93.6 nm for the PIM-1/ZIF-8 composite membrane in the z-direction. Atoms within 1 nm from both ends of the membranes were tethered to their initial z-coordinates to prevent their drifting due to the created concentration gradient. To create the initial configurations of gas molecules, they were randomly placed into the void spaces on both sides of the membranes. Single-component H2 and CH4 as well as H2/CH4 mixture simulations were performed. In all simulations, the concentration of the gas molecules in the inlet control region (feed) was maintained at their experimentally measured molecular density at 5 bar and 300 K,[53] which were 0.1203 and 0.1217 molecules/nm3 for H2 and CH4, respectively. Thus, the considered mixture corresponds to almost an equimolar feed composition. Outlet gas concentration was set to vacuum. Periodic boundary conditions were applied in all directions. Simulations were run in the NVT ensemble and the temperature of the systems was fixed at 300 K using a Nosé–Hoover thermostat. The thermostat coupling constant was set to 0.1 ps. Separate thermostats were used for fluid molecules and individual membrane components to prevent asymmetric thermalization in the simulation box due to hot solute–cold solvent effect.[54] ZIF-8 and PIM-1 were modeled with all-atom and united atom flexible force fields, respectively. Lennard–Jones (LJ) parameters and partial charges for the ZIF-8 and PIM-1 atoms as well as details of intramolecular force field terms can be found elsewhere.[29] CH4 and H2 were modeled with united atom force fields, both implementing an uncharged single LJ site, with parameters taken from Martin et al.[55] and Frost et al.,[56] respectively. Particle Mesh Ewald method was employed to account for long-range electrostatics interactions. A 1.2 nm cutoff distance was used for the LJ and the real part of the Ewald sum. LJ cross-term parameters for the interactions between membrane atoms and gas molecules were tuned to capture the magnitude of the available experimental H2/CH4 permselectivity data in literature for the individual PIM-1 and ZIF-8 membranes (Table ).[37,57−65] These refined parameters given in Table S2 were also used to study the PIM-1/ZIF-8 composite. The equations of motion were integrated with a 1 fs time step using a Verlet scheme. Single-component permeation simulations of H2 and CH4 through the PIM-1 and ZIF-8 membranes and the PIM-1/ZIF-8 composite membrane were run for 1 μs each. In addition, a H2/CH4 mixture separation simulation through the PIM-1/ZIF-8 membrane was also run for 1 μs. Single-component and mixture simulation results were reported for the last 200 ns. We should emphasize that normally about 100 ns of simulation is sufficient to achieve a steady-state diffusion; however, the runs were extended to the microsecond scale with the aim of demonstrating the computational feasibility of the CGD-MD simulations without suffering any feed depletion issues as discussed elsewhere.[33]

Table 1

Comparison of Experimental and Simulated Ideal H2/CH4 Permselectivities in PIM-1 and ZIF-8 Membranesab

PIM-1	10.4[58]	6.47[64]	5.25[59]	5.42[57]	8.37[65]	6.27ba
ZIF-8	13.0[61]	4.61[62]	12.5[60]	4.86[37]	4.63[63]	5.12ba

Superscripts refer to the references that the experimental data are taken from.

Based on single-component CGD-MD simulations in this work.

Superscripts refer to the references that tn class="Chemical">he experimental data are taken from.

The flux of n class="Chemical">H2 and CH4 gases along the z-direction (J) was calculated by counting the net number of molecules that cross an xy-plane located at the center of the membrane and dividing it by the simulation time (t) and the cross-sectional area of the membrane (A)where N+ and N– are the number of H2 or CH4 molecules that cross the xy-plane in the +z-direction (i.e., feed to permeate) and the −z-direction (i.e., permeate to feed), respectively. The fluxes were then used to calculate gas permeabilities (Table S3) and H2/CH4 permselectivities.

Residence time probability distributions of the H2 and CH4 molecules within the membranes were obtained by calculating the time spent by individual molecules within 1 nm long bins along the z-direction. One important clarification that needs to be made here is that the residence time of a molecule in a bin was determined regardless of the direction (i.e., positive and negative z-directions) it has entered and left the bin. The residence time probability distributions were then used to calculate the mean residence times in each bin.

Results and Discussion

The CGD-MD approach consists of creating a concentration gradient across the membrane that facilitates the transport of gases. As such, it is essential that the concentrations of molecules are maintained at their target values at the inlet and outlet of the membrane. Figure S2 shows the concentration of H2 and CH4 molecules in the inlet control and outlet control regions of the membrane systems studied. In all systems simulated, the CGD-MD method succeeds in keeping the concentration of the gases very close to the target values.

Table reports the comparison between tn class="Chemical">he simulated ideal H2/CH4 permselectivities in PIM-1 and ZIF-8 membranes and the available experimental ideal permselectivities. The ideal permselectivity is calculated by taking the ratio of single-component permeabilities. The experimental permselectivites for both PIM-1 and ZIF-8 vary within a broad range; however, the simulated ideal permselectivities lie within this range and are in good agreement with several of the reported permselectivity data. Overall, both membranes are H2 selective, since H2 permeability is greater than that for CH4 in both PIM-1 and ZIF-8.

Figure shows the density profiles of H2 and CH4 molecules in the PIM-1 and ZIF-8 membranes along the direction of gas flow (i.e., z-direction). The H2 density values are close to each other in both membranes, that is, H2 is adsorbed in similar amounts in PIM-1 and ZIF-8. On the other hand, CH4 adsorption is higher in PIM-1 compared to that in ZIF-8. Furthermore, the density of CH4 is higher than that of H2 in both membranes, indicating a stronger adsorption of CH4 compared to that of H2 in PIM-1 and ZIF-8. As a result, while the H2 density decreases almost linearly due to the concentration gradient, there is a sharp increase in the CH4 density at the entrance of the membranes before gradually decreasing. The H2 molecules quickly permeate through the membranes and do not exhibit a density increase at the entrance of the membrane compared to its density in the feed. Conversely, a strongly adsorbed CH4 exhibits a much higher density compared to its density in the feed.

Figure 2

Density profiles of the single-component H2 and CH4 along the (a) PIM-1 and (b) ZIF-8 membranes. Dashed lines correspond to the location of membrane surfaces. Snapshots from single-component CH4 and H2 CGD-MD simulations are given to visually guide concentration gradients along the PIM-1 (c and e, respectively) and ZIF-8 (d and f, respectively) membranes. Green and yellow spheres represent CH4 and H2 molecules, respectively.

Figure shows the residence time analyses of the single-component H2 and CH4 permeation along the PIM-1 and ZIF-8 membranes. In the PIM-1 membrane, H2 exhibits a relatively flat profile, whereas in the ZIF-8 membrane, H2 molecules exhibit longer residence times, thanks to the presence of cages. In contrast to that of H2, the residence time profile of CH4 in PIM-1 exhibits large variations due to the relatively stronger adsorption of CH4 in PIM-1. Local structural fluctuations in the PIM-1 membrane lead to a rather inhomogeneous residence time profile along the membrane. For instance, the relatively high mean residence time for CH4 at the z-coordinate of around 23 nm for PIM-1 indicates the presence of a relatively large cavity around this location, where CH4 molecules spend more time. In the ZIF-8 membrane, both H2 and CH4 mean residence time profiles exhibit an up-and-down pattern, reflecting ZIF-8’s repeated structure composed of large cages connected by narrow apertures; i.e., mean residence times of the molecules are longer in the cages but shorter near the apertures. In both membranes, CH4 residence times are longer than H2 residence times.

Figure 3

Mean residence times of H2 and CH4 along the z-direction of the (a) PIM-1 and (b) ZIF-8 membranes. Each point corresponds to a mean residence time within a 1 nm wide bin. Dashed lines correspond to the location of membrane surfaces.

Mean residence times of H2 and n class="Chemical">CH4 along the z-direction of the (a) PIM-1 and (b) ZIF-8 membranes. Each point corresponds to a mean residence time within a 1 nm wide bin. Dashed lines correspond to the location of membrane surfaces.

After simulating single-component H2 and CH4 permeations in PIM-1 and ZIF8 membranes, we investigated their single-component and mixture permeations in the composite PIM-1/ZIF-8 membrane. Figure compares the resulting density profiles of H2 and CH4. The H2 density almost linearly decreases along the composite membrane in both cases. The reason for such a linearity along the entire composite structure, despite H2 permeating through different structures, is that H2 is adsorbed in similar quantities in PIM-1 and ZIF-8 (Figure a,b). In contrast, the CH4 density decreases along the first PIM-1 slab, then exhibits a drop at the PIM-1/ZIF-8 interface, and continues to decrease along ZIF-8 before showing a sharp increase at the interface between ZIF-8 and the second PIM-1 slab. The drop in CH4 density after the first PIM-1 slab is a consequence of the fact that ZIF-8 adsorbs less CH4 compared to PIM-1, and the jump after the onset of the second PIM-1 slab is because PIM-1 adsorbs more CH4 compared to ZIF-8 (Figure a,b). There are also some quantitative differences in the z-density profiles of H2 and CH4 for the single-component and mixture simulations. It can be seen that H2 density is lower throughout the entire membrane in the mixture simulation when compared to the single-component case. This is because CH4 molecules occupy some of the spaces that were previously available only for the H2 molecules in the single-component simulation. On the other hand, the difference between the density profiles of CH4 for the single-component and mixture simulations is less significant.

Figure 4

Density profiles of (a) CH4 and (b) H2 along the z-direction of the composite PIM-1/ZIF-8 membrane in single-component and mixture simulations. Dashed lines correspond to the location of membrane surfaces. Snapshots from single-component (c) CH4 and (d) H2 CGD-MD simulations are given to visually guide concentration gradients along the composite PIM-1/ZIF-8 membrane. Green and yellow spheres represent CH4 and H2 molecules, respectively.

Table compares the computed n class="Chemical">H2/CH4 permselectivities in the composite PIM-1/ZIF-8 membrane. The H2/CH4 permselectivity obtained from the CGD-MD simulation for the mixture in the composite PIM-1/ZIF-8 membrane is lower than the ideal H2/CH4 permselectivity calculated based on single-component permeabilities. The mixture value is deemed more accurate because the simulation of the mixture takes into account the interactions between the H2 and CH4 molecules. As previously explained, in the mixture simulation CH4 molecules displace H2 molecules in the membrane (Figure ), and as a consequence, H2 permeability decreases in the mixture simulation while the permeability of CH4 in the single-component and mixture simulations are almost the same (Table S2). The deviation between the ideal and mixture H2/CH4 permselectivities clearly demonstrates the importance of taking interactions between two different gas species into account when predicting permselectivities in membranes, which is overlooked in the calculation of the ideal permselectivity.

Table 2

Comparison of H2/CH4 Permselectivities in the Composite PIM-1/ZIF-8 Membrane Obtained by Different Methods

H2/CH4 permselectivity	method used
4.61	ideal permselectivity based on single-component CGD-MD simulations.
3.36	permselectivity from mixture CGD-MD simulation.
5.51	ideal permselectivity based on serial resistance model[28] using CGD-MD-computed H2 and CH4 permselectivities in PIM-1 and ZIF-8 (Table 1).

While a discrepancy between the ideal and mixture H2/CH4 permselectivities is expected to a certain degree, the CGD-MD simulations of the composite PIM-1/ZIF-8 membrane reveal a much more critical issue on the effect of the polymer/MOF interface in predicting the permselectivity of the composite membrane. In principle, the ideal permselectivity of an MMM is expected to lie between the permselectivities of its constituent materials. That is, the ideal H2/CH4 permselectivity of the PIM-1/ZIF-8 membrane should be between the ideal H2/CH4 permselectivities of PIM-1 and ZIF-8. One of the macroscopic permeation models that we can conveniently use to predict the permselectivity of the composite PIM-1/ZIF-8 membrane, based on H2 and CH4 permeabilities in PIM-1 and ZIF-8 (Table S2), is the serial resistance model,[28] which is defined aswhere Peff is the effective permeability of the composite for a given gas and and ϕ are the individual permeability and volume fraction of the constituents of the composite membrane. Indeed, when the serial resistance model is employed, the predicted ideal H2/CH4 permselectivity of the composite PIM-1/ZIF-8 membrane (Table ) is between those for PIM-1 and ZIF-8 (Table ). However, the permselectivity from the serial mode is about 20% larger than the ideal H2/CH4 permselectivity calculated based on single-component permeabilities from the CGD-MD simulations of the composite PIM-1/ZIF-8 membrane, i.e., 5.51 vs 4.61, respectively. This difference can be attributed to the nonselective microvoids that exist at the PIM-1/ZIF-8 interfaces, resulting in interfacial resistances. In this case, the interfacial resistance between PIM-1 and ZIF-8 is relatively significant, such that the CGD-MD-computed ideal H2/CH4 permselectivity of the composite PIM-1/ZIF-8 membrane is even lower than that obtained for ZIF-8 (Table ). While the CGD-MD simulations can directly include the effect of such interfacial resistances in predicting the permselectivity of a composite material for a given gas separation, the serial resistance model does not take such effects into account.

Conclusions

In this study, we report microsecond-long CGD-MD simulations of H2 and CH4 transport in PIM-1 and ZIF-8 membranes as well as their permselectivity in the composite PIM-1/ZIF-8 membrane. The CGD-MD method allowed us to map the density and mean residence time profiles of the H2 and CH4 gases along the membranes while these two gases diffuse under a concentration gradient. Furthermore, we directly computed the flux of the permeating gases through the membranes from the CGD-MD simulations and evidenced the effect of the interfaces on the H2/CH4 permselectivity in the composite PIM-1/ZIF-8 membrane. We found that the presence of nonselective void spaces between PIM-1 and ZIF-8 in the composite PIM-1/ZIF-8 membrane induces a decrease of the H2/CH4 permselectivity by about 20% compared to the ideal permselectivity estimated by the application of macroscopic model on the data obtained individually for ZIF-8 and PIM-1. The CGD-MD simulations carried out with an accurate description of the polymer/MOF interfaces allow the determination of the magnitude of such a deviation, thus paving the way for a more critical use of macroscopic models to predict the performances of MMMs. This work provides a first unambiguous proof that interfaces play a crucial role in the gas transport mechanism in polymer/MOF composites. This makes questionable the extensive use of macroscopic models to predict the performances of MMMs, since these models do not take into account either defects at the interface or interactions between guest molecules in gas mixtures, factors that have a big impact in the performance of the MMM.