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 mem-branes 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 H 2 / CO 2 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 H 2 / CH 4 and CO 2 /CH 4 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 CO 2 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 H 2 and CH 4 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 H 2 /CH 4 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 m 2 /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, H 2 /CH 4 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 H 2 and CH 4 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
2.1. 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 zdirection (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  1.
2.2. 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. Singlecomponent H 2 and CH 4 as well as H 2 /CH 4 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/nm 3 for H 2 and CH 4 , 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 Nose−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 CH 4 and H 2 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 H 2 /CH 4 permselectivity data in literature for the individual PIM-1 and ZIF-8 membranes (Table 1). 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 H 2 and CH 4 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 H 2 /CH 4 mixture separation simulation through the PIM-1/ZIF-8 membrane was also run for 1 μs. Singlecomponent 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 The flux of H 2 and CH 4 gases along the z-direction (J z ) 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 xy ) where N i + and N i − are the number of H 2 or CH 4 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 H 2 /CH 4 permselectivities.
Residence time probability distributions of the H 2 and CH 4 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 Chemistry of Materials pubs.acs.org/cm Article 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 H 2 and CH 4 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 1 reports the comparison between the simulated ideal H 2 /CH 4 permselectivities in PIM-1 and ZIF-8 membranes and the available experimental ideal permselectivities. The ideal permselectivity is calculated by taking the ratio of singlecomponent 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 H 2 selective,  Chemistry of Materials pubs.acs.org/cm Article since H 2 permeability is greater than that for CH 4 in both PIM-1 and ZIF-8. Figure 2 shows the density profiles of H 2 and CH 4 molecules in the PIM-1 and ZIF-8 membranes along the direction of gas flow (i.e., z-direction). The H 2 density values are close to each other in both membranes, that is, H 2 is adsorbed in similar amounts in PIM-1 and ZIF-8. On the other hand, CH 4 adsorption is higher in PIM-1 compared to that in ZIF-8. Furthermore, the density of CH 4 is higher than that of H 2 in both membranes, indicating a stronger adsorption of CH 4 compared to that of H 2 in PIM-1 and ZIF-8. As a result, while the H 2 density decreases almost linearly due to the concentration gradient, there is a sharp increase in the CH 4 density at the entrance of the membranes before gradually decreasing. The H 2 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 CH 4 exhibits a much higher density compared to its density in the feed. Figure 3 shows the residence time analyses of the singlecomponent H 2 and CH 4 permeation along the PIM-1 and ZIF-8 membranes. In the PIM-1 membrane, H 2 exhibits a relatively flat profile, whereas in the ZIF-8 membrane, H 2 molecules exhibit longer residence times, thanks to the presence of cages. In contrast to that of H 2 , the residence time profile of CH 4 in PIM-1 exhibits large variations due to the relatively stronger adsorption of CH 4 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 CH 4 at the zcoordinate of around 23 nm for PIM-1 indicates the presence of a relatively large cavity around this location, where CH 4 molecules spend more time. In the ZIF-8 membrane, both H 2 and CH 4 mean residence time profiles exhibit an up-and-down pattern, reflecting ZIF-8's repeated structure composed of large Chemistry of Materials pubs.acs.org/cm Article 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, CH 4 residence times are longer than H 2 residence times. After simulating single-component H 2 and CH 4 permeations in PIM-1 and ZIF8 membranes, we investigated their singlecomponent and mixture permeations in the composite PIM-1/ ZIF-8 membrane. Figure 4 compares the resulting density profiles of H 2 and CH 4 . The H 2 density almost linearly decreases along the composite membrane in both cases. The reason for such a linearity along the entire composite structure, despite H 2 permeating through different structures, is that H 2 is adsorbed in similar quantities in PIM-1 and ZIF-8 ( Figure  2a,b). In contrast, the CH 4 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 CH 4 density after the first PIM-1 slab is a consequence of the fact that ZIF-8 adsorbs less CH 4 compared to PIM-1, and the jump after the onset of the second PIM-1 slab is because PIM-1 adsorbs more CH 4 compared to ZIF-8 (Figure 2a,b). There are also some quantitative differences in the z-density profiles of H 2 and CH 4 for the single-component and mixture simulations. It can be seen that H 2 density is lower throughout the entire membrane in the mixture simulation when compared to the single-component case. This is because CH 4 molecules occupy some of the spaces that were previously available only for the H 2 molecules in the single-component simulation. On the other hand, the difference between the density profiles of CH 4 for the single-component and mixture simulations is less significant. Table 2 compares the computed H 2 /CH 4 permselectivities in the composite PIM-1/ZIF-8 membrane. The H 2 /CH 4 permselectivity obtained from the CGD-MD simulation for the mixture in the composite PIM-1/ZIF-8 membrane is lower than the ideal H 2 /CH 4 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 H 2 and CH 4 molecules. As previously explained, in the mixture simulation CH 4 molecules displace H 2 molecules in the membrane ( Figure  4), and as a consequence, H 2 permeability decreases in the mixture simulation while the permeability of CH 4 in the singlecomponent and mixture simulations are almost the same (Table S2). The deviation between the ideal and mixture H 2 / CH 4 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.
While a discrepancy between the ideal and mixture H 2 /CH 4 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 H 2 /CH 4 permselectivity of the PIM-1/ZIF-8 membrane should be between the ideal H 2 /CH 4 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 H 2 and CH 4 permeabilities in PIM-1 and ZIF-8 (Table S2), is the serial resistance model, 28 which is defined as where P eff is the effective permeability of the composite for a given gas and Pi and ϕ i 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 H 2 /CH 4 permselectivity of the composite PIM-1/ZIF-8 membrane (Table 2) is between those for PIM-1 and ZIF-8 (Table 1). However, the permselectivity from the serial mode is about 20% larger than the ideal H 2 /CH 4 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 H 2 /CH 4 permselectivity of the composite PIM-1/ZIF-8 membrane is even lower than that obtained for ZIF-8 (Table 1). 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 H 2 and CH 4 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 H 2 and CH 4 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 H 2 /CH 4 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 H 2 /CH 4 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 ideal permselectivity based on serial resistance model 28 using CGD-MD-computed H 2 and CH 4 permselectivities in PIM-1 and ZIF-8 ( Table 1).
Chemistry of Materials pubs.acs.org/cm Article 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.
Schematic representation of CGD-MD simulation setup; CGD-MD-specific parameters used in simulations; variation of CH 4 and H 2 concentrations as a function of simulation time in the control regions; tuned Lennard−Jones cross parameters; and permeabilities of H 2 and CH 4 (PDF)