Adsorption space for microporous polymers with diverse adsorbate species

The enormous number of combinations of adsorbing molecules and porous materials that exist is known as adsorption space. The adsorption space for microporous polymers has not yet been systematically explored, especially when compared with efforts for crystalline adsorbents. We report molecular simulation data for the adsorptive and structural properties of polymers of intrinsic microporosity with a diverse set of adsorbate species with 345 distinct adsorption isotherms and over 240,000 fresh and swollen structures. These structures and isotherms were obtained using a sorption-relaxation technique that accounts for the critical role of flexibility of the polymeric adsorbents. This enables us to introduce a set of correlations that can estimate adsorbent swelling and fractional free volume dilation as a function of adsorbate uptake based on readily characterized properties. The separation selectivity of the 276 distinct binary molecular pairs in our data is reported and high-performing adsorbent systems are identified.


INTRODUCTION
In one form or another, chemical separations have an unavoidable presence in the functioning of nearly every sector of the chemical process industry. By some estimates 10-15% of global energy consumption stems from separation processes 1 , which has yielded opportunities for emergent low-energy adsorption-based technologies that could result in an annual savings of billions of US dollars in the United States 2,3 . One avenue for progress toward this challenge is the replacement of traditional energy-intensive industrial processes such as distillation with lower energy separations driven by microporous adsorbent materials. Alongside the demand for separation applications, microporous adsorbents can have potential impact in other areas, such as chemical storage 4 or catalysis 5 .
To address the wide-reaching nature of applications mentioned above various families of microporous adsorbent materials have been developed. Metal-organic frameworks 6 , covalent-organic frameworks 7 , and zeolites 8 are common types of adsorbents with crystalline order, whereas hyper-cross-linked polymers 9 , conjugated microporous polymers 10 , carbonaceous materials 11 , and polymers of intrinsic microporosity (PIMs) 12 are examples of adsorbents with amorphous structures. A common feature to all of these families is that thousands or even millions of structures can be envisioned on the basis of chemical modifications or systematic combinations of constituent components. However, substantial experimental exploration of a library of these structures tends to be impractical, especially in cases where materials synthesis is complex. This challenge can potentially be addressed through molecular simulation screening studies before the development of synthetic chemistry methods 13,14 , which can also generate the data needed for machine learning approaches 15 .
PIMs are defined by contorted and ladder-like backbones that inefficiently pack and produce high-free volume microporous polymeric adsorbents. Through creative design of monomeric units many remarkable PIMs have been produced. For example, Rose et al. 16 presented 2-dimensional PIM-TMN-Trip membranes that surpass the performance of poly(1-(trimethylsilyl)-1-propyne) (PTMSP), a polymeric material currently used commercially. Recent work by Thompson et al. 17 demonstrated the use of PIM structures for the fractionation of light shale-based crude oil. An attractive characteristic of PIMs is their potential to be scaled up to meet the demands of large volume industrial separations. They can be solution processed, readily modified by a range of organic chemical reactions, and present in different geometries (e.g., films or hollow fibers). Despite their potential, information that addresses a large and diverse collection of PIM-adsorbate pairs is not currently available. One explanation for the lack of these studies are the difficulties associated with modeling amorphous microporous polymers. The simulation techniques needed to construct models of amorphous polymers have been discussed by Kupgan et al. 18 . The demands of these techniques mean that there is no exhaustive structure or property database of amorphous microporous polymers that rivals the size of the databases of crystalline materials 19,20 . An initial step toward providing access to existing structures was reported recently by Thyagarajan and Sholl, who made available a database of 205 atomically detailed amorphous materials collected from previous studies 21 .
It is not only the structure of amorphous polymers that complicate modeling of adsorption in PIMs. Microporous polymers of all kinds are influenced by varying degrees of adsorbateinduced framework rearrangement. This limits property predictions from traditional molecular simulation methods because the modeling of adsorption including polymer flexibility is a nontrivial task 22 . Depending on the pressure range, adsorbates, and adsorbent species studied, employing a simulation method that accounts for polymer dynamics, swelling, and plasticization can be essential for making quantitative predictions. To highlight a specific example, an adsorbent's separation selectivity can be severely reduced by swelling/plasticization, an effect that is entirely unaccounted for with a traditional rigid framework approximation 23,24 .
Here, we provide an extensive dataset of amorphous polymeric adsorbents under different loading conditions (over 240,000 structural conformations), including 345 unary adsorption isotherms accounting for polymer flexibility, and an evaluation of the binary selectivity of 4140 PIM-mixture combinations. We evaluate the effect that systematic chemical modification of the adsorbent monomeric units has on adsorption loading and show how several adsorptive properties and structural feature rearrangements can be determined from simple relationships with reasonable accuracy.

RESULTS
Adsorption with polymer flexibility Molecular simulations that account for framework flexibility for adsorption predictions can quantify structural features that may be otherwise difficult to obtain experimentally, e.g., determining a polymer's fractional free volume (FFV) at elevated adsorbate pressures. FFV can be determined experimentally using positron annihilation lifetime spectroscopy, but to the best of our knowledge this technique has never been applied with an adsorbed phase at significantly swollen conditions. FFV can also be probed experimentally with X-ray scattering, inverse gas chromatography, and N 2 adsorption 25 ; however, demonstrations of extensibility for swollen adsorbents are limited. In addition to FFV, simulations can measure the expansion of geometric surface area (SA geo ) and induced volumetric swelling (V/V o ) associated with adsorption. To provide context for the results to follow, the behavior of carboxyPIM-1 upon adsorption of propene, acetaldehyde, butane, hexane, and dipropyl ether is discussed. Figure 1a shows five adsorption isotherms in carboxyPIM-1 at 300 K. In general, the larger molar mass adsorbates have stronger adsorption in the Henry's Law regime than the lighter species. As relative pressure increases, however, the heavier adsorbates display saturation-like behavior, whereas the loadings of the lighter species continue to increase without clear saturation. Figure 1b provides swelling (V/V o )-pressure curves that correspond to the state points of Fig. 1a. We report swelling as a relative value with reference to the adsorbate-free polymer rather than a percent change; a relative swelling of 1.1 is equivalent to a 10% volume change. Figure 1c, d shows that the framework flexibility predicted in the highlighted adsorbate-carboxyPIM-1 combinations is substantial, with values ranging from 0.208 to 0.508 and 362 to 2221 m 2 g −1 for FFV dilation and SA geo expansion, respectively. The full set of loading, swelling, FFV dilation, and surface area expansion data for every adsorbate-adsorbent pair other than propan-1-ol is provided in the repository listed under Data availability.
Although MC/MD simulations are effective for determining a series of state points on an adsorption isotherm, it is useful to fit the simulated data with an isotherm model. We used the dualmode sorption (DMS) model [26][27][28] to fit the 345 adsorption isotherms obtained from our simulations. The accuracy of the DMS model fit is shown in Fig. 2a. The DMS functional form gives a mean absolute deviation (MAD) of 0.16 mmol g −1 . As a demonstration, the adsorption isotherms from Fig. 1 have been replotted with their DMS fits in Fig. 2b. The fitted DMS model parameters are reported in the provided database. Validation for the sorptionrelaxation procedure is presented in Supplementary Fig. 1 for propane and propene. The discrepancy between traditional Grand Canonical Monte Carlo (GCMC) assuming a rigid polymer structure and the sorption-relaxation methods for n-hexane with carboxyPIM-1 is shown in Supplementary Fig. 2, which highlights the considerable additional uptake from swelling.

Correlations for sorption-induced rearrangement
To progress towards adsorbent design criteria for polymers it is essential to have methods to quantify polymer restructuring. Experimentally, measurements of macroscopic swelling are possible with spectroscopy 29,30 or in situ ellipsometry 31 , although these approaches are not yet widespread in the development of high-free volume polymeric membrane materials. Our simulations provide a set of over 240,000 polymer conformations exhibiting a range of adsorption-induced swelling. This large set of data allows us to consider relationships that can be used to estimate polymeric framework swelling from the adsorption isotherm and the connection between swelling and FFV dilation. We found that the V/V o -Q curve is convex with two regimes 24 . The initial swelling regime, V/V o < 1.15, follows a quadratic relationship described by where K VQ is a fitted constant specific to the adsorbate-adsorbent system. For V/V o > 1.15 swelling is approximately linear, following with γ VQ also being a system specific fitted constant. Figure 3a illustrates that this piecewise pair of equations is effective in representing the data for all adsorbing species. The MAD is 0.004 across all samples, which is appreciably less than the 95% confidence intervals associated with the data points. The propensity to swell the polymer framework is generally correlated with the molar mass of the adsorbate species. Figure 3b shows an example of this observation. For a given adsorbent, the values of K VQ and γ VQ can be estimated by and with T C being the adsorbate's critical temperature and MM the adsorbate's molar mass. We find that the leading coefficients of Eqs. (3) and (4) can be related to adsorbent structural features by and ξ ads ¼ 4:57ðSA geo V He Þ À1=4 (6) where SA geo is the geometric surface area (m 2 g −1 ) and V He the helium void volume (cm 3 g −1 ) of the adsorbate-free polymer. Physically, the product of SA geo and V He is a descriptor of the available adsorbent pore space in the empty polymer. It is useful to note that these equations could be used in both computational and experimental studies because they rely only on quantities accessible by routine characterization of the fresh polymer. In Fig. 5a the FFV dilation is plotted as a function of relative swelling. At low-swelling values (i.e., where elastic deformation is expected to occur) the FFV dilation is occurring to a degree approximately three times the volumetric swelling. This allows for the modification of the Eq. (3) (restricted to V/V o < 1.15) to represent the relative FFV dilation as Dilation of this nature results in a slight shift of the pore size distribution (PSD) and a peak intensity change dictated by the expanded FFV. This type of PSD for V/V o = 1.11 is shown in Fig. 5b. In contrast, for a material that is 11% swollen relative to the static framework the simulated results display a mean peak position at a larger value with greater intensity and broadened distribution. The PSD in the 33% swollen (severely plasticized) material shows the same behavior to a greater extent. This discrepancy reveals that significant free volume restructuring is occurring in a nonaffine nature.

Evaluating binary separation performance
To report separation performance of PIMs for adsorption of binary mixtures, we exploit the fact that the binary adsorption selectivity (α) in the dilute limit is exactly equal to the ratio of the singlecomponent Henry's constants 32 . At non-dilute conditions the restructuring of polymeric materials generally leads to a reduction of α 33 . The Henry's constants used for calculating α were determined using the Widom insertion method 34 and their distribution is reported in Supplementary Fig. 9. Additionally, a box and whisker plot of the dilute heats of adsorption are provided in Supplementary Fig. 10. We characterize adsorbent selectivity performance for a given binary mixture as the maximum of α and 1/α, where a selective adsorbent displays values greater than 1. The spread across all PIMs is reported in Fig. 6a, where the best α for each molecular pair is shown above the diagonal and the worst α below the diagonal. This is by far the largest set of data for mixture adsorption in PIMs that is currently available.
To expand the adsorbent selectivity analysis, a simple relationship for predicting the typical α for arbitrary binary pairs, similar to the functional form reported by Tang et al. 35,36 that was based on a relationship described by Freeman 37 , was explored. Specifically, we assessed the applicability of predicting the dilute limit selectivity using where T C;A À T C;B is the difference of adsorbate critical temperatures. By fitting the simulated median α values for each adsorbate pair we obtain a scaling factor specific to PIMs: N pim = 0.043 K −1 . This equation fits our simulated data well, as shown in Fig. 7, although the variations associated with the full range of polymer architectures are considerable.

DISCUSSION
We have reported simulated adsorption data for several high-free volume PIMs in combination with 24 diverse adsorbate species. Importantly, the simulated results are obtained with a sorptionrelaxation technique that overcomes the 'rigid' framework approximation, which allows for quantification of adsorbateinduced polymeric adsorbent restructuring. The results of these simulations include 345 adsorption, swelling, FFV dilation, and surface area expansion isotherms, and a database of over 240,000 fresh and swollen polymer structures. These results provide insight for the role that species diversity plays in adsorption. However, the observations are contingent upon the quality of the modeling approach. For the set of PIMs studied in this work, structure factors from wide-angle X-ray scattering, adsorptive properties, and mechanical behavior have been previously validated 16,[38][39][40] . This is essential to the physical relevance of the results produced.
Ensuring that all of these properties are accurately represented is key to applying the sorption-relaxation approach. As an example, if the PIM models display mechanical characteristics that are too weak because of an incorrect interatomic potential, then the sorption-induced swelling will be overpredicted. Further information regarding modeling details and implications can be found in References 16,38-40 as well as Supplementary Discussion 1. In Fig. 1a the degree of swelling is 1.0-1.05 at low relative pressures for all species. A similar result was found in 2018 by Kupgan et al. 23 for PIM-1, which led to the conclusion that this regime could be simulated with either GCMC/MD or traditional rigid MC with approximately equivalent accuracy for the lighter species of their work (methane, carbon dioxide, nitrogen, oxygen, and hydrogen). However, Fig. 2 shows that at a pressure of P/P o = 0.1 a n-hexane loading of~2.0 mmol g −1 results in framework dilation of~1.13; in other words, beyond the applicability of the rigid framework assumption. The challenge of utilizing PIMs at such conditions would be that the deliverable capacity of a plasticized material severely suffers with each adsorption-desorption cycle, i.e., depreciating hysteresis behavior. An attractive future approach is to focus on developing PIMs that can be regenerated quickly and at low cost, to take advantage of the large uptake due to their propensity to swell. Figure 2c, d shows that the degree of relative change of adsorbent structural features follows the order of volumetric swelling < FFV dilation < SA geo expansion. This has implications for tuning adsorbent performance since FFV is known to play a dominant role in adsorbate mass transport 41 and SA geo has a considerable influence on adsorption thermodynamics (e.g., enthalpy of adsorption). The observed behaviors are widely consistent across all adsorbate-adsorbent pairs assessed in this work. Propan-1-ol was omitted from the analysis in Fig. 2 because the majority of state points were found to excessively swell or did not reach convergence even after hundreds of MC/MD iterations. This is likely a result of the sigmoidal isotherm shape that this species displays in PIMs 42 and/or the molecular clustering that can occur for small primary alcohols 43 . An empirical DMS model was found to fit the adsorption isotherm data with a MAD of 0.16 mmol g −1 . Polymer restructuring is not accounted for in the DMS functional form and therefore it is somewhat surprising that this choice is so effective in fitting our diverse dataset. We used the data from our molecular simulations to develop a set of correlations for PIMs that give good estimates for swelling and FFV dilation at moderate swelling (V/V o < 1.15) and reasonable predictions across the full range of swelling. These correlations rely only on readily accessible properties of the polymer, SA geo and V He . SA geo can be easily calculated computationally if a structural model is available using open-source software (e.g., Poreblazer 44 or zeo++ 45 ) or by subtracting 250 m 2 g −1 from the experimentally measured Brunauer-Emmett-Teller surface area (SA BET ) as shown by Hart and Colina 46 . Pore volume is measurable using helium pycnometry or can be computationally calculated using the same software mentioned above. To connect the quality of our reported correlations with the experimental sensitivity of swelling measurements, we highlight a recent experimental work by Ogieglo et al. 47 that showed a variability of 2-4% in swelling induced by CO 2 in PIM-1 (with V/V o < 1.1) based on the initial polymer film thickness. This variability is similar to the spread of the initial swelling data points in Fig. 4. We emphasize that the method introduced relies only on characterization of the adsorbent in the unswollen fresh state without adsorbate present (SA geo and/or V He ). These adsorbent properties are related to separation thermodynamics and mass transport, which can be tied to the well-known permeability-selectivity tradeoff that polymers exhibit 48,49 . It seems likely that these relationships could be relevant to other  Fig. 6 Evaluation of room temperature separation (α) performance for 15 PIMs with 253 adsorbate binary pairs. a α matrix (log 10 scale) with the best performing separations across the 15 PIMs shown above the diagonal and the worst performance shown below. The bluer/ redder the color the better the separation performance. b Histogram of the frequency a given PIM's separation performance is above the average of the 15 PIMs. c Histogram of the frequency a given PIM appears as the top performer for a binary adsorbate mixture. families of high-free volume glassy polymers, or even more generally to high-free volume organic glasses.
We found that adsorbent restructuring upon adsorption is significantly nonuniform (nonaffine), which manifests in the relative change in FFV dilation occurring at a greater degree than the overall adsorbent swelling. This occurs because during swelling free volume elements are dilating in relation to the adjusted polymer chain positions but at the same time additional free volume elements and pore connectivity are being created. It is this second effect that drives the change of FFV to be greater than that of the swelling behavior. Regardless of the approach for measuring volumetric swelling (whether it is through the above correlations or direct experimental measurement) the results of Fig. 5a should be equally extendable.
We also quantified the dilute separation performance of 4140 adsorbent-binary mixture pairs, allowing us to identify highperforming adsorbents. By comparing symmetric entries in Fig. 6a across the diagonal we see that for the majority of binary mixtures a range of α spanning orders of magnitude is available among the PIMs. Challenging separations of adsorbate pairs, e.g., with similar sizes and/or chemistry are also identified. soPIM-1 was observed to be the overall top performer across the set of PIMs with the greatest number of best separations capabilities and being consistently above average. This appears to be related to the adjacent position of the sulfonyl group to the nitrile, which creates strongly interacting electrostatic sites within the polymer that enhances separation capabilities for mixtures containing polar species 50 . We also found that estimates for typical selectivity behavior for PIMs follows a simple relationship based on the adsorbates' critical temperature. It should be noted that there are several instances where the preferential species is split among the PIMs studied, in other words, some of the frameworks favor species A over B and the rest B over A. These instances can be seen in Supplementary Fig. 11. There are also seven gas mixtures where all of the PIMs assessed prefer that opposite species to the one predicted with Eq. (8). The variation in selectivity is relatively large and this correlation should be viewed as a qualitative baseline to compare against adsorbent-specific selectivities. The field of microporous polymeric materials has progressed tremendously since the introduction of the original correlation, i.e., 20 years ago and before the synthesis of the first PIM. The usefulness of this relationship lies on the possibility of predicting the typical PIM selectivity and not the selectivity of a specific PIM material. Moreover, the dilute limit separation capabilities are known to have poor extension to loading conditions where swelling and plasticization occur. An example from Kupgan et al. showed a 20% reduction in CO 2 /CH 4 selectivity when PIM-1 is simulated at plasticizing conditions 23 . In general terms, this reduction is associated with regions of expanded volume in a structure that has plasticized and/or swelled. A significant complication for making quantitative predictions about multicomponent selectivity in these systems at non-dilute loadings is that the degree of swelling must be a function of the activity of all of the adsorbing species. An interesting challenge for future work will be to explore whether information on the swelling from multiple singlecomponent adsorption isotherms can be combined to predict that full range of behaviors that arise from adsorbed mixtures. Additionally, a highly swollen pore structure can provide spatial accessibility to large species (e.g., n-octane) and result in regions of competitive adsorption that were previously size excluded.
Finally, to assess whether privileged adsorbents 36 exist among the examples we considered two metrics as reported in Fig. 6b, c, namely the frequency at which a given adsorbent performs better than the average α and the frequency with which a given adsorbent gave the best α from the set. A similar analysis with an alternative metric that combines selectivity and gas uptake is demonstrated in Supplementary Fig. 12. While the notion of a privileged adsorbent is not defined by specific quantitative criteria, it is certainly appropriate for soPIM-1 in the set of the adsorbents studied. This particular polymer appears above average for 74% of the separations considered and is the top performer among all the PIMs we considered for 68% of them. PIM-1 displays low comparative performance (23% above average and 0% best separations). Even if soPIM-1 and soPIM-0 are not considered in the analysis of Fig. 6c, PIM-1 still exhibits a similar comparative performance (23% above average and 0% best separations). This finding reiterates that there is significant membrane design space available through chemical modifications that can be explored in parallel with efforts to develop PIMs composed of more complex monomers 25,51,52 .
This work provides an expanded understanding of the interplay between adsorbent diversity, adsorbate diversity, and polymer restructuring as a function of adsorbate loading for microporous polymeric adsorbents. Although computationally laborious, our molecular simulations have provided insights that would otherwise be challenging to obtain experimentally. We provided a database of over 240,000 fresh and swollen PIM structures and 345 distinct adsorption isotherms. The single-component isotherms and adsorbent restructuring and selectivity relationships produced will be valuable for materials design and predicting adsorption behavior in the very large adsorption space associated with polymeric adsorbents. The initial success of the correlations reported here for single-component adsorption suggests that similar relationships may ultimately be useful for describing the more challenging situation of multicomponent adsorption in flexible polymers.

Polymer structure preparation and modeling parameters
All polymer structures assessed in this work were previously built 53,54 using the open-source Polymatic software 55 . Briefly, Polymatic simulates 'computational polymerization' by starting with a low-density simulation cell packed with a defined number of repeating units and creating bonds using distance and orientation criteria. Between bond searches, MD simulations are performed to allow the growing polymer chains to sample the simulation cell space and find other unreacted species. After a specified number of bonds have been added the polymer systems are further equilibrated to production quality samples using methods based on the work of Theodorou and Suter 56 , which was later adapted by Larsen et al. for PIMs 38 . Essentially, a set of NPT and NVT ensemble MD simulations are performed at increasing temperatures/pressures, followed by incremental cooling/depressurization, and ultimately equilibration at standard temperature and pressure (300 K and 1 atm). Conditions for each simulation step are provided in Supplementary Table 1. Producing physically representative polymer configurations is vital to making meaningful adsorption predictions 18 . The monomer units of the 15 polymer chemistries are shown in Fig. 8 with the associated naming convention.
The interatomic potential applied to the polymeric adsorbents was the General Amber Force Field 57 for bonded terms and the united-atom variation of the Transferable Potential for Phase Equilibria (TraPPE-UA) for nonbonded terms [58][59][60][61][62] . Electrostatic partial point charges were derived using the restrained electrostatic potential fitting method 63 . This combination of force fields, in conjunction with the structure generation procedure referenced above, has been well-validated for modeling a range of PIM materials by reproducing wide-angle X-ray scattering spectra, adsorptive properties, and tensile moduli 16,[38][39][40] . Force field information for each monomer unit illustrated in Fig. 8 is openly available in the Colina group's force field database. Porous features of the equilibrated polymer samples studied in this work were characterized with Poreblazer (v3.0.2) 44 and are reported in Supplementary Table 2. The adsorbates of interest were modeled as fully flexible species with TraPPE/TraPPE-UA (see Supplementary Table 3 for more details). Characteristics of each adsorbent polymer conformation are summarized in Supplementary Tables 4-6. The TraPPE force field is based on accurate bulk phase properties of a large set of molecules and is appropriate for our adsorbate species of interest. Longrange electrostatics were calculated using Ewald summation 64 for MC simulations and the particle-particle particle-mesh (p3m) 65

Molecular simulation details
Adsorption simulations were performed at 300 K with a combined GCMC and MD technique; sometimes referred to as the sorption-relaxation approach. This method has been previously applied to study sorptioninduced polymer rearrangement in PIMs 23,50,66 and is appropriate for evaluating adsorption performance under swelling/plasticizing conditions. It is shown above that an approach that incorporates structural rearrangement (such as sorption-relaxation) is essential for the PIMs we consider at relatively high loadings; there are many state points evaluated that would be significantly misrepresented, generally as an underprediction of adsorption loading, without accounting for adsorbent restructuring. This paper reports the utilization of a computational technique that accounts for framework flexibility for a study of appreciable size, which required >300 µs of MD simulations and 3 × 10 9 MC cycles (roughly equivalent to 6 months of wall time with 2500 CPUs on the University of Florida's HiPerGator 2.0 system).
Adsorption-relaxation is carried out in multiple iterations where GCMC simulations are used to predict adsorbate loading and MD simulations are employed to account for system rearrangement. The simulated adsorption isotherms were produced up to the vapor pressure of each adsorbate species (Supplementary Table 7), unless excessive swelling is observed; in other words, samples displaying uncontrolled dilation (continuous volume expansion that tears apart the polymer chains). This excessive swelling was generally observed only for species with a high solubility at a P/P o of 0.5-1.0. Each GCMC portion of the sorption-relaxation approach was performed using RASPA 67 for 10 5 cycles, where a cycle is the number of MC steps equal to the greater value of 20 and the number of adsorbate species present. MD simulations were carried out with the LAMMPS software package (version 3Mar2020) 68 for 1 ns with barostat conditions corresponding to the adsorbate's bulk phase pressure. For adsorbates modeled with united atoms the 1 ns MD simulation was split in two parts for the purpose of efficiency; the first 100 ps used a 1 fs timestep and the remaining 900 ps used a 2 fs timestep. For adsorbates that were modeled with explicit hydrogens (thiol and alcohol containing molecules) we used a 1 fs timestep for the full MD trajectory. An assumption in this sorptionrelaxation approach is related to adsorbate kinetic accessibility. In traditional MC simulations, such as the GCMC simulations we used, the entire simulation cell is considered sampleable and thus all pores are assumed to be accessible. It is possible in some polymer conformations, however, that a portion of pores may in reality be kinetically isolated. We elected to not block any pores to account for this effect during the sorption-relaxation simulations.

Statistical analysis
Polymer models were built with five distinct replicates, and thus all data reported are averaged over the independent simulations of each system with error bars denoting a 95% confidence interval. After every GCMC-MD iteration the system is checked for convergence based on the observed averages for volume, energy, and adsorbate loading. The number of iterations to reach convergence varied among the examples we considered, but on average a single pressure point required~25 iterations (i.e., 25 ns of MD and 2.5 × 10 6 GCMC cycles). The number of GCMC-MD iterations needed for convergence at a defined state point is not known a priori. To address this important issue in an automated way the number of GCMC/MD iterations is controlled by adopting MC convergence statistics as described previously 50,66 based on the work of Gelman and Rubin 69 . A Gelman-Rubin factor (R) is calculated using by considering statistics of the last m GCMC/MD iterations (m = 5 in this work) that contain a number of data points n. The defining feature of Eq. (9) is the last ratio that compares the mean variance of each iteration (B) to the variance of the means between the iterations (W). Conceptually, R approaches 1 for a completely converged system (i.e., simulations of infinite length). In practice R is larger than 1 and decreases as the system approaches convergence. We continued our simulations until R < 1.5. In test calculations, this threshold meant that symmetric fluctuations around average loading values were observed and increasing simulation duration would mainly serve to narrow confidence intervals.