A high-throughput screening of metal – organic framework based membranes for biogas upgrading †

Applications of biomethane as a source of renewable energy and transport fuel rely heavily on successful implementation of puri ﬁ cation methods capable of removing undesirable impurities from biogas and increasing its calori ﬁ c content. Metal – organic frameworks (MOFs) are competitive candidates for biogas upgrading due to a versatile range of attractive physical and chemical properties which can be utilised in membrane materials. In this work, we present a high-throughput computational screening methodology for e ﬃ cient identi ﬁ cation of MOF structures with promising gas separation performance. The proposed screening strategy is based on initial structural analysis and predictions of the single-component permeation of CO 2 , CH 4 and H 2 S from adsorption and di ﬀ usion calculations at in ﬁ nite dilution. The identi ﬁ ed top performing candidates are subject to further analysis of their gas separation performance at the operating conditions of 10 bar and 298 K, using grand canonical Monte Carlo and equilibrium molecular dynamics simulations on equimolar CO 2 /CH 4 and H 2 S/CH 4 mixtures. The Henry constant for the adsorption of H 2 O was also calculated to determine the hydrophobicity of MOF structures, as the presence of H 2 O often leads to membrane instability and performance limitations. For the considered gas mixtures, the top MOF candidates exhibit superior separation capabilities over polymer-, zeolite-, and mixed matrix-based membranes as indicated by the predicted values of selectivity and permeability. The proposed screening protocol o ﬀ ers a powerful tool for the rational design of novel MOFs for biogas upgrading. available from the literature. 35 In this work, a diverse set of 7909 MOF structures are found to be suitable for high-throughput screening. The RASPA so  ware package 46 is employed to calculate the self-di ﬀ usivity ( D 0 ), Henry constant ( K 0 ) and heat of adsorption ( Q 0st ) of CO 2 , H 2 S, CH 4 and H 2 O at in  nite dilution. K 0 is determined using the Widom insertion method at 298 K, using a total of 100 000 cycles. 47 For calculations of D 0 , EMD simulations in the NVT ensemble are employed. A time step of 1 fs is used, and the number of equilibration and production cycles are set to 10 000 and 100 000, respectively. Heating is controlled using a Nos´e – Hoover thermostat. 47 To mimic in  nite dilution, adsorbate – adsorbate interactions are neglected and only framework – adsorbate interactions are considered. The statistical accuracy of in  nite dilution conditions is improved by using 30 gas molecules per simulation. The initial positions of the gas molecules are determined using 1000 Monte Carlo initialization cycles. The mean square displacement (MSD) of the gas molecules is computed and the self-di ﬀ usivity is obtained using the gradient of the MSD via Einstein ’ s relation. 47 Based on previous work, we exclude materials where D 0 < 10 (cid:1) 8 cm 2 s (cid:1) 1 as these values approach the limit of accurate description of molecular di ﬀ usion by EMD simulations. 48 The gas


Introduction
Agricultural, industrial and municipal wastes are some of the many feedstocks which, in the absence of O 2 , can be broken down to form biogasa mixture of gases composed of CH 4 (50-65%), CO 2 (35-50%) and small amounts of trace gases such as H 2 S, N 2 , H 2 and NH 3 . 1 The composition of biogas is largely dependent on the waste substrate and the conditions under which substrate digestion occurs. 2 Raw biogas has a high caloric content between 5000-7000 kcal m À3 , making it ideal for heat and electricity generation. 3 Through separation of its components, biogas can be upgraded to provide an enriched source of CH 4 , further improving its caloric content for utilisation in the transport sector. This process is essential as high levels of CO 2 lower the caloric content of the fuel. In addition, small concentrations of H 2 S (>50 ppm) have been shown to corrode metal, making raw biogas unsuitable for standard combustion engines. 4 In industry, there are many techniques that can be employed to upgrade biogas including absorption, adsorption and cryogenic distillation methods. Absorption methods involve the capture of impurities by dissolving them in a solvent at elevated pressure, or by reacting them with an inorganic/amine based solution. 5 Adsorption methods make use of molecular sieves, activated carbon and zeolites to separate out biogas due to differences in affinity between the biogas components and the adsorbent bed. The contaminants are usually adsorbed at high pressure from the biogas stream, leaving CH 4 to pass over the adsorbent bed. The bed can then be regenerated by reducing the pressure which desorbs the contaminants. These are the core principles underlying pressureswing adsorption methods used to separate gaseous mixtures. 6 Similar principles can also be employed using changes in other stimuli such as temperature. Finally, as each biogas component condenses at different temperatures and pressures, one can use distillation methods to separate contaminants. 7 These methods can provide high CH 4 purity typically between 90-99%, however, they can be quite energy intensive, suffering from high investment, operation and/or maintenance costs. 8 One of the most promising and cheaper alternatives is to use membrane materials to separate gas components through differences in permeation. Membranes can be synthesised from various materials, including organic polymers, molecular sieves, palladium alloys, amorphous silicas, ceramics and zeolites. 9 Polymer based membranes currently dominate the membrane market as they are cheap, easily scalable, can be fabricated into hollow bers or sheets and have been more extensively studied. However, they suffer from short lifetime, low thermal and chemical stability, and low selectivity. 10 Furthermore, the work of Robeson revealed a trade-off relationship between selectivity and permeability in polymer membranes which limits their overall separation performance. 11,12 Inorganic based membranes tend to produce higher selectivity and possess greater thermal and chemical stability than polymer membranes making them suitable for separations that require harsh conditions. 13 However, they typically have higher production costs and can have reproducibility issues when attempting to create uniform, defect-free membranes. 14 While a lot of progress has been made in the development of synthesis methods for existing membranes, [15][16][17] the discovery of new materials which can rival the separation performance of existing membranes is also a key area of research. [18][19][20] One example of these materials are metal-organic frameworks (MOFs): a rapidly growing hybrid class of porous materials that are selfassembled from inorganic clusters (known as secondary building units) and organic ligands. Due to the diverse nature of their building blocks, MOFs can offer a versatile range of attractive physical and chemical properties such as high porosity and large surface area. 21 Furthermore, these properties are highly separation performance and stability of MOF membranes. 44 In this work, we apply a multi-stage screening methodology to focus on CO 2 /CH 4 and H 2 S/CH 4 separations across a range of hydrophobic MOF membranes from the CSD database. Top selected candidates encompass high selectivity and permeability for CO 2 and H 2 S which is essential for the biogas upgrading process.
Additionally, a set of membranes with improved moisture stability due to their enhanced hydrophobic character has been identied. The initial ltering of porous candidates is based on the differences in their geometric properties. Next, a time efficient estimation of the ideal CO 2 /CH 4 and H 2 S/CH 4 selectivity and permeability is carried out at innite dilution conditions. This is achieved by calculating the Henry coefficients using a grand canonical Monte Carlo (GCMC) approach and the diffusion coefficients extracted from equilibrium molecular dynamics (EMD) calculations. We also explore how different input parameters, used to calculate partial charges via the extended charge equilibration (EQeq) scheme, 45 can affect the physical nature and the overall separation properties of MOFs. The top hydrophobic MOFs found above the Robeson upper bound are then subjected to calculations of the performance prediction under working conditions. We perform equimolar CO 2 /CH 4 and H 2 S/CH 4 GCMC and EMD simulations to predict mixture selectivity and permeability at 298 K and 10 bar. The best selected membranes are compared to literature zeolites, MOFs and MMMs to assess their overall performance and viability for biogas upgrading. Finally, we examine the structural properties of 1183 hydrophobic membranes that lie above the Robeson limit. This enables us to draw conclusions about the relationships which lead to high separation performance, providing useful insight into the design of future membranes.

Computational details
A set of geometric criteria, detailed in the ESI, † is applied to the CSD MOF data set (67 675 structures) to remove structures with non-accessible surface area and pore diameters smaller than the individual components of biogas. Solvent molecules are removed from each structure using a python script readily available from the literature. 35 In this work, a diverse set of 7909 MOF structures are found to be suitable for high-throughput screening. The RASPA soware package 46 is employed to calculate the self-diffusivity (D 0 ), Henry constant (K 0 ) and heat of adsorption (Q 0 st ) of CO 2 , H 2 S, CH 4 and H 2 O at innite dilution. K 0 is determined using the Widom insertion method at 298 K, using a total of 100 000 cycles. 47 For calculations of D 0 , EMD simulations in the NVT ensemble are employed. A time step of 1 fs is used, and the number of equilibration and production cycles are set to 10 000 and 100 000, respectively. Heating is controlled using a Nosé-Hoover thermostat. 47 To mimic innite dilution, adsorbate-adsorbate interactions are neglected and only framework-adsorbate interactions are considered. The statistical accuracy of innite dilution conditions is improved by using 30 gas molecules per simulation. The initial positions of the gas molecules are determined using 1000 Monte Carlo initialization cycles. The mean square displacement (MSD) of the gas molecules is computed and the self-diffusivity is obtained using the gradient of the MSD via Einstein's relation. 47 Based on previous work, we exclude materials where D 0 < 10 À8 cm 2 s À1 as these values approach the limit of accurate description of molecular diffusion by EMD simulations. 48 The gas permeability for each component i at innite dilution (P 0 i ) is calculated from the product of self-diffusivity and the Henry constant, i.e. P 0 i ¼ D 0 i Â K 0 i . The ratios K 0 i / K 0 j , D 0 i /D 0 j and P 0 i /P 0 j of two components are then used to determine the ideal adsorption (S 0 ads,i/j ), diffusion (S 0 diff,i/j ) and membrane selectivity (S 0 mem,i/ j ), respectively. A hydrophobic subset of MOFs is identied by comparing K 0 H2O of each membrane against the hydrophobic ZIF-8 structure (K 0 H2O ¼ 5.00 Â 10 À6 mol kg À1 Pa À1 ). A membrane is considered hydrophobic if the value of K 0 H2O is smaller or equal to that of ZIF-8. This approach has been used previously to identify hydrophobic MOFs for hydrocarbon separations and toxic industrial chemical capture. 49,50 Only hydrophobic membranes above the 2008 Robeson bound 12 that are selective to both H 2 S and CO 2 and possess P 0 CO2 and P 0 H2S greater than 2.5 Â 10 7 barrer are considered in the nal phase of simulations where their performance at working conditions of 10 bar and 298 K is analysed. The value of 2.5 Â 10 7 barrer is chosen to reduce the total number of structures and does not reect a particular target permeability. For GCMC simulations, the number of initialization and production cycles is set to 10 000 each (20 000 total). Gas phase fugacities are calculated using the Peng-Robinson equation of state. 51 Due to the high cost of the EMD simulations, the LAMMPS soware package is used to enable parallelisation of the calculations. 52 The initial states of each EMD simulation are determined from the loadings predicted in the GCMC simulations. In contrast to calculations at innite dilution, both hostadsorbate and adsorbate-adsorbate interactions are included. A time step of 1 fs is employed with simulations running for a total of 1 ns (50% equilibration). Atomic coordinates for each adsorbate are recorded every 5 fs and are used to calculate the MSD and self-diffusivity via Einstein's relation. A minimum of 5 trajectories are used to compute the mixture self-diffusivity of each component. The gas permeability in the binary mixture is re-evaluated using and f i are the gas permeability, uptake, self-diffusion coefficient and fugacity of component i, respectively, and F is the fractional pore volume of the membrane. 53 The membrane selectivity (S mix mem,i/ j ) can then be calculated from the ratio of the gas permeabilities for components i and j of the binary mixture as P mix i /P mix j .

Force-eld description and validation
The force-elds used to describe intermolecular interactions are formed from a Lennard-Jones (LJ) potential with a cut-off of 12.8 A. Long-range coulombic interactions are handled using the Ewald summation technique. 54 Cross-terms are calculated using the Lorentz-Berthelot mixing rules. 55,56 Frameworks are assumed to be rigid with atoms xed in their crystallographic positions and modelled using the Universal Force-Field (UFF). 57 Atomic partial charges are assigned to framework atoms using the EQeq scheme available within the RASPA soware package. 45 A list of oxidation states used in EQeq calculations can be found in the ESI. † CO 2 and CH 4 molecules are represented by rigid 3-site and single united-atom models, respectively, using LJ parameters and partial charges from the TraPPE force-eld. 58,59 H 2 O is represented by a rigid 4-site model using the TIP4P-Ew model 60 and H 2 S is represented by a rigid 4-site model using parameters derived by Kristóf and co-workers. 61 A full list of force-eld parameters may be found in Tables S1 and S2. † For MOFs, so called "off-the-shelf" force-elds such as UFF are commonly employed to study adsorption and diffusion properties. These are typically combined with force-elds designed specically for adsorbates to reproduce thermophysical properties, such as the TraPPE and TIP force-elds. This computational setup has been used extensively over the past 15 years and has been shown to give good agreement across a wide range of adsorbates and MOFs. 62-64 Furthermore, they have been previously validated for reproducing experimental membrane properties in a variety of CO 2 and CH 4 based separations. 65,66 In addition to these previous studies, we have also performed similar calculations which show good agreement between experimental and simulated values for permeability of CO 2 and CH 4 ( Fig. 1top).
Studies of H 2 S adsorption and separation are less documented than for other small molecules such as CO 2 and CH 4 . 69 In many MOFs, the formation of strong metal-sulphur bonds can lead to poor regeneration of the empty host or even displacement of linkers causing loss of permanent porosity. Despite these issues, a small number of MOFs from the MIL family display H 2 S adsorption isotherms over large pressure ranges, which are suitable for validation of the H 2 S model. Our results on MIL-47(V) indicate good agreement with experimental results (Fig. 1 bottom) and similar agreement to simulations using a 3-site H 2 S model employed by Hamon and co-workers. 67 Good agreement between theory and experiment is also achieved for MIL-68(Al). 68 An evaluation test for H 2 S permeability over a greater range of MOF membranes is not feasible with the current availability of experimental data. A perfect, rigid, solvent-and defect-free MOF crystal structure is assumed in all calculations, thus providing an upper limit of permeation properties achievable in these materials. Whilst this approximation has been used ubiquitously in the literature, some membranes may be sensitive to guest activation leading to structural changes. Additional chemical and mechanical stability tests of the candidates identied in a high-throughput screening may be also required.
One limitation of the force-elds used in this work is that strong interactions between adsorbates and open metal sites are not captured, leading to underprediction of the separation performance in these materials. Studies that have attempted to parametrise force-elds using quantum chemical calculations have successfully shown that this method can correctly capture host-guest interactions around open metal sites. 70,71 While this approach has been used at the highthroughput scale previously, it is a very costly process which can only target a small subset of MOFs possessing similar structural motifs (e.g. metal paddlewheels). 72,73 Finally, a number of materials, e.g. some MOFs containing openmetal sites, are omitted from this screening due to poor separation performance which may arise from limitations in the employed models. As our aim is to target a range of MOFs with a diverse set of structural properties, general force-elds such as UFF will give the best compromise of accuracy and cost.

Accuracy in evaluation of partial charges on metal sites
In the early stages of the screening, the distribution of partial charges were evaluated for 7909 MOFs using the EQeq method. Unlike the original Qeq method reported by Rappe et al., 74 this method uses a neutral oxidation state Taylor expansion for ligand atoms, and a higher oxidation state Taylor expansion for metals. In some systems, the choice of Taylor expansion can be the difference between convergence to a physical set of partial charges or an unphysical set. In previous work, Ongari and co-workers compared the partial charges obtained from different charge equilibration methods against high quality density functional theory-derived DDEC charges for a large set of MOFs. 75 They showed that smaller deviations between the Qeq methods and DDEC charges were obtained when using a higher oxidation state for the metal centres present in the MOFs. Furthermore, using the data from DDEC charge distributions, they showed that by imposing an upper and lower atomic charge limit of +3 and À2, respectively, fewer structures containing unphysical charges were obtained from using the higher oxidation state Taylor expansion. In this work, MOFs with physical partial charges are obtained by applying two sets of criteria. Firstly, MOFs containing atoms with a partial charge greater than +3.5 or less than À2 were removed. This criterion is closely based on the condition proposed by Ongari and co-workers, however, we nd that a slight buffer added to the upper limit is useful for structures with highly charged metal centres containing Zr(IV). We nd that this criterion alone is sufficient for identifying metals with unphysical charges but less so for ligand atoms such as carbon or hydrogen ( Fig. S1 and S2 of the ESI †). Therefore, we apply a stricter charge criteria on C and H atoms to only allow MOFs with C and H charges between +1 and À1. A total of 800 MOFs are removed from the screening leaving 7109 MOFs to be taken forward into the rst phase of highthroughput screening.
To demonstrate how the choice of Taylor expansion may inuence the results, Fig. 2 compares Henry coefficients for CO 2 and H 2 S (K 0 CO2 and K 0 H2S ) in MOFs for which a converged EQeq partial charge evaluation was achieved using both neutral and higher oxidation states of the metal atoms. For both CO 2 and H 2 S, using the higher oxidation state Taylor expansion (x-axis) provides a narrower and more chemically meaningful spread of values for Henry coefficients. This is because the neutral oxidation state expansion has a higher propensity to provide structures with unphysical charges. 75 When using the neutral oxidation state Taylor expansion, we analyse the charge distributions in 76 MOFs which possess a K 0 CO2 . 10 3 mol kg À1 Pa À1 , (i.e. any value greater than the largest K 0 CO2 predicted from using the higher oxidation state expansion). We nd that 50% of these MOFs contain Li, Na or K sites with charges ranging from +1.34 to +6.75, and 30% contain Nd sites with charges ranging from +2.89 to +46.4. Consequently, the unphysical Coulomb interactions between CO 2 and these metals centres lead to large overestimation of K 0 CO2 . Using the higher oxidation state expansion, the Li, Na and K atoms in these MOFs have charges ranging from +0.55 to +0.71. Additionally, the Nd charges ranging between +1.45 and +2.46.
A total of 4418 MOFs (68%) have a larger K 0 CO2 when using the higher oxidation state expansion instead of the neutral oxidation state expansion (Fig. 2 le). However, in 81% of these MOFs, this increase in K 0 CO2 is relatively small (1 to 10 times larger). Although this change is small, it is worth noting that the adsorption selectivity obtained from the ratio of the CO 2 and CH 4 Henry constants can be directly inuenced by switching from one Taylor expansion to another. Similar Fig. 2 Comparison of the Henry coefficients for CO 2 (red) and H 2 S (green) in MOFs, predicted by using different sets of partial charges on metal sites. The y-axis corresponds to partial charges calculated using a neutral oxidation state on metal atoms, while the xaxis corresponds to partial charges calculated using a higher oxidation state on metal atoms. A full list of oxidation states may be found in the ESI. † trends are also observed for H 2 S with 4013 MOFs (62%) possessing a larger K 0 H2S when using the higher oxidation state expansion. In 88% of these MOFs, K 0 H2S is 1 to 10 times larger compared to that calculated using the neutral oxidation state expansion. In general, these results demonstrate the importance of using a higher oxidation state Taylor expansion for metal atoms for generating a more accurate and physical set of partial charges.

Permeability at innite dilution
Within the solution-diffusion model, the process of gas permeation across a membrane can be broken down into three steps: adsorption onto the membrane surface, diffusion through the membrane pore network, and desorption from the membrane surface. Therefore, it is expected that differences in host-guest and guestguest interactions for different biogas components will strongly inuence the rate at which a component permeates through the membrane. For an initial insight, calculations at very low concentrations (in the Henry region) are useful for comparing differences in permeation as these are solely inuenced by the strength of host-guest interactions. Fig. 3 compares the diffusion coefficients and permeability of 6768 MOF membranes for CH 4 , CO 2 and H 2 S at innite dilution conditions. We note to the reader that from the initial 7109 structures, 341 membranes are excluded as they possess diffusion coefficients smaller than 10 À8 cm 2 s À1 .
For CH 4 gas, there appears to be a positive correlation in which structures possessing large P 0 CH4 are likely to possess large K 0 CH4 and D 0 CH4 (Fig. 3a). For 5264 membranes, D 0 CH4 lies between 10 À3 and 10 À4 cm 2 s À1 , resulting in fast rates of diffusion which we attribute to weak CH 4 -MOF interactions. Although K 0 CH4 values range between 8.91 Â 10 À3 mol kg À1 Pa À1 and 1.58 Â 10 À7 mol kg À1 Pa À1 , 81% of MOFs possess K 0 CH4 \3:00 Â 10 À5 mol kg À1 Pa À1 also indicating weak CH 4 -MOF interactions. Unlike CH 4 gas, both CO 2 and H 2 S interact with MOF membranes through a combination of van der Waals and electrostatic interactions. Therefore, on average, we expect CO 2 -MOF and H 2 S-MOF interactions to be stronger than CH 4 -MOF interactions. We see evidence of this in the range of K 0 CO2 and K 0 H2S values which extend further than K 0 CH4 to values as large as 10 3 mol kg À1 Pa À1 (Fig. 3b and c). Furthermore, 76% of MOFs possess a K 0 CO2 between 10 À3 mol kg À1 Pa À1 and 10 À5 mol kg À1 Pa À1 (70% for K 0 H2S ). Therefore, on average these gases diffuse slower than CH 4 through each MOF membrane. While 21% of MOFs possess D 0 CH4 \2:00 Â 10 À4 cm 2 s À1 , this increases to 63% and 66% when considering D 0 CO2 and D 0 H2S . The combination of these factors lead to a more negative looking correlation as presented in Fig. 3b and    separation performance of these membranes, Robeson plots were constructed, as shown in Fig. 4. For CO 2 /CH 4 , the 1991 and 2008 Robeson upper bounds 11,12 are included to illustrate the performance of these MOFs relative to some of the best available organic polymer membranes. 1499 MOFs are found to be more selective for CO 2 gas with S 0 mem;CO2=CH4 ranging between 1 and 26.9 and P 0 CO2 values ranging between 4.10 Â 10 4 barrer and 1.56 Â 10 8 barrer. As seen in this study, MOFs have the potential to overcome the Robeson upper bound limit corresponding to the selectivity-permeability trade-off problem found in organic polymer membranes. Highlighted in red and yellow in Fig. 4 are a total of 1183 MOFs that exceed the 2008 Robeson limit for CO 2 /CH 4 separation. For H 2 S/CH 4 separation, only 14 MOFs were found to be more selective for CH 4 than H 2 S. The remaining 1671 possess S 0 mem;H2S=CH4 and P 0 H2S ranging from 1.03 to 82.3 and from 7.23 Â 10 4 to 1.98 Â 10 9 barrer, respectively. Due to the limited availability of H 2 S permeability data, an upper limit for the H 2 S/CH 4 separation is not established. Therefore, the structures highlighted in red and yellow in Fig. 4b are the same as those shown in Fig. 4a and are highlighted for comparative purposes.
Due to the large computational cost associated with modelling molecular dynamics at high pressure, we focus our study on MOF candidates that show the most promise at innite dilution conditions. Hence, we consider MOF structures above the upper bound that are selective to CO 2 and H 2 S with P 0 CO2 and P 0 H2S greater than 2.50 Â 10 7 barrer, leaving a total of 42 structures. Seven of these structures (CSD reference codes: CAFSUY, CAFTAF, COMDOY, CUMDIY, LOD-POL, OHAKEO and SETQAL) are reported to collapse in response to activation or possess degrees of structural exibility when exposed to different stimuli. [76][77][78][79][80][81] Therefore, the rigid lattice approximation used in our calculations will not provide an accurate representation of their structure used under working conditions. These MOFs were removed from the screening leaving a total of 35 unique MOFs represented by yellow circles in Fig. 4.

Separations at working conditions
For the promising membrane candidates, we employ GCMC and EMD simulations to study equimolar binary mixture separations of CO 2 /CH 4 and H 2 S/CH 4 at 10 bar and 298 K. Unlike calculations at innite dilution, these simulations include the role of adsorbate-adsorbate interactions in biogas separation at operating conditions. One MOF structure (REFCODE: HONCIY) was removed from the screening as an accurate estimation of D mix CO2 and D mix CH4 was not possible within the simulation timescale.
A large decrease in gas permeability P mix CH4 , P mix CO2 and P mix H2S is observed at working conditions, with values ranging between 10 1 and 10 5 barrer (Fig. 5). For 22 of the 34 MOFs, the values of both S mix mem;CO2=CH4 and S mix mem;H2S=CH4 are greater than their ideal selectivity at innite dilution. With the inclusion of guest-guest interactions and increased gas loading under working conditions, the rate of diffusion for each component in the gas mixture is reduced leading to slower permeation rates. Furthermore, each component now competes against one another for adsorption sites which is not captured in the ideal calculations at innite dilution. Under working conditions, the separation performance of 20 MOFs deteriorates such that they now lie below the 1991 and 2008 Robeson bounds (Fig. 5top). Of the remaining 14 MOFs, 6 lie between the 1991 and 2008 bounds and 8 lie above both bounds. These materials possess S mix mem;CO2=CH4 between 3.09 and 12.5, P mix CO2 between 2.87 Â 10 4 barrer and 1.24 Â 10 5 barrer, S mix mem;H2S=CH4 between 25.0 and 123 and P mix H2S between 2.72 Â 10 3 barrer and 1.03 Â 10 5 barrer.

Top candidates
In order to rank the performance of each MOF based on the values of P mix CO2 , P mix H2S , S mix CO2=CH4 , and S mix H2S=CH4 , a scoring system has to be generated by weighting each of these parameters. This is because much larger concentrations of CO 2 need to be removed from a biogas stream than H 2 S, which may affect the rankings of one MOF to another. However, the list of top structures is strongly correlated with the type and size of the applied weighting constraints, which in some cases identied top structures that fall below the Robeson upper bounds. Furthermore, it is possible that a weighting system of such design may introduce bias into the selection of the top candidates. One possible solution to this problem is to convert each property into a descriptor that is passed through a machine-learning algorithm. This may generate rankings in a more unbiased manner based on the descriptors and problem alone, however this lies outside the scope of this publication. For the remainder of this discussion, the top 8 MOFs are referred to as the top candidates identied by our high-throughput screening and correspond to those which exceed the 2008 CO 2 /CH 4 Robeson limit (Fig. 6). While this does not give a quantitative rank for each material, it does highlight a set of MOF membranes that tend to perform better than other candidates in terms of permeability and selectivity. Void analysis of the top candidates indicates that they all possess very similar structural characteristics such as one-dimensional pore channels (see ESI †). The geometric and physical properties summarised in Table 1 indicate that in these materials, the largest cavity diameter (LCD) is typically smaller than 6 A, accessible surface area (ASA) ranges from 164 to 1018 m 2 g À1 , and void fraction (VF) is typically less than 0.5. Of the top 8 candidates, the MOF with CSD reference code QOKCID has the smallest LCD resulting in the lowest gravimetric and volumetric storage capacities for CO 2 and H 2 S. In contrast, the RUVBER MOF exhibits the largest LCD with capacity for more CO 2 and H 2 S storage than most top candidates. Pores with large diameters typically possess a greater proportion of sites with weak host-guest interactions, affecting their selectivity for contaminants. Despite this, RUVBER is still able to selectively  Diffusion coefficients D mix CO2 , D mix CH4 and D mix H2S of the top candidates range from 4.38 Â 10 À5 cm 2 s À1 to 7.87 Â 10 À7 cm 2 s À1 . The average diffusion selectivity S mix diff for each binary mixture is fairly close to unity, ranging between 0.53 and 1.14 for CO 2 /CH 4 mixtures and between 0.60 and 3.25 for H 2 S/CH 4 mixtures. In the CO 2 /CH 4 mixtures, CO 2 molecules interact more strongly with the framework than CH 4 resulting in a reduced degree of freedom that is generally re-ected in their smaller self-diffusion coefficients. In the H 2 S/CH 4 mixtures, the large gas loadings and high S mix ads;H2S=CH4 result in each CH 4 molecule being surrounded by a large number of contaminant molecules. This strongly inuences the diffusion rate of CH 4 such that both H 2 S and CH 4 molecules possess similar rates of diffusion through the pore channels. In general, D mix CO2 $ D mix H2S because each MOF possesses a lower S mix ads;CO2=CH4 value than S mix ads;H2S=CH4 . Consequently, a larger concentration of CH 4 is found in the pores during the CO 2 /CH 4 MD simulations which can collide with nearby CO 2 molecules to help increase their diffusion rate.
We next compare the top MOF candidates to those of other porous membranes reported in the literature, which include pure zeolite-, MOF-and ZIF-based membranes, in addition to a number of novel MMMs ( Table 2). Note that many studies report CO 2 separations over a range of conditions which may not be exactly the same as those used in biogas upgrading. However, these studies are still useful for evaluating the separation performance of the top MOF candidates identied in the proposed screening. Some of the largest values of S mix mem were reported by Carreon 84 and Shi 85 for SAPO-34 membranes, and many of the reported MOF-and ZIF-based membranes possess P mix CO2 between 10 4 barrer and 10 5 barrer and S mix mem ranging between 1 and 25. For many MOFs, S mix mem appears to improve in composite MMMs at the cost of lowering the value of P mix CO2 which typically ranges from 10s to 1000s of barrer. However, it is clear from the studies reported in Table 2 that changes to the polymer and MOF loading can signicantly alter the permeability, resulting in an additional level of control over the permeation properties. In comparison, it appears that the top MOFs membranes identied in this work generally possess larger S mix mem and P mix CO2 values. Furthermore, as these materials lie above the Robeson limit, their separation performance is competitive with traditional polymer membranes. As previously mentioned, our simulations assume a perfect, defect-free crystal structure, in which the values of P mix CO2 and S mix mem represent an upper performance limit. Despite this, the top MOFs in this work appear to show promise in separating biogas contaminants. Further insight which assesses the top membranes as ller particles in MMMs would be benecial to establish whether their separation performance can be enhanced further. However, this consideration goes beyond the scope of the current study.

Structure-performance relationships
One of the main advantages of working with large sets of structures is the ability to investigate the relationships between their performance and physical properties, which are crucial in the rational design of future materials. Generally, the separation performance for a class of materials cannot be linked to an individual structural, physical or chemical property; it requires understanding of a complex interplay between all of them. Fig. 4 shows that 1183 MOFs can surpass the selectivity-permeability trade-off relationship described by Robeson's upper bound. By comparing the structural properties of the MOFs considered in this study, we aim to answer three closely linked questions: (1) are there any common properties shared by these 1183 structures? (2) How do these properties account for the observed separation performance? And (3) are there any relationships between the structural properties and the separation performance that will allow us to construct design principles for future MOFs?
Histogram distributions shown in Fig. 7 for different structural properties of the 1183 MOFs indicate that many of these contain small pore sizes, with 80% of the structures having a pore-limiting diameter (PLD) less than 6 A and 81% possessing LCD below 8 A. The ultra-microporous nature of many of these materials also leads to small ASAs with 68% found to possess an ASA less than 1250 m 2 g À1 . Of the remaining membranes that possess an ASA greater than 1250 m 2 g À1 , only 42% possess an LCD greater than 8 A. Large pores with high surface area can make it more difficult to separate mixtures through differences in binding affinity, making them undesirable for these types of separation. The porosity distribution in Fig. 7a indicates that MOFs with low VF appear more frequently, with 42% having a VF < 0.4 and 61% < 0.5. Analysis of the metals in each structure indicates that 29% of the structures are constructed from zincbased nodes, forming the most common metal ahead of copper (14%), cadmium (14%) and cobalt (11%). Furthermore, these four metals were exclusively found in 7 of the 10 top MOFs. In a recent study by Altintas and co-workers, it was shown that these four metals also frequently appear in the top membranes predicted for CH 4 /H 2 separations. 65 Finally, the density distribution in Fig. 7f shows that 64% of MOFs above the upper bound have densities between 1 g cm À3 and 1.5 g cm À3 . In summary, our calculations at innite dilution suggest that MOFs with a PLD between 3.8 A and 6 A, LCD between 4 A and 8 A, ASA less than 1250 m 2 g À1 , VF between 0.3 and 0.5 and density between 1 g cm À3 and 1.5 g cm À3 are likely to appear above the Robeson upper bound and exhibit desirable separation performance.
It is useful to compare these results with the structural property distributions from low-performing materials that lie below the 2008 Robeson bound. If structures below the upper bound also have similar distributions, then this suggests that these properties are not of signicance and do not correlate with high performance. Interestingly, it is found that only 10% of membranes below the upper bound possess similar structural properties which all lie in the same ranges that were summarised at the end of the previous paragraph. Two of the most signicant differences appear to be the PLD and LCD. In the low-performing subset of membranes, 42% possess a PLD smaller than 6 A, and 43% have a LCD smaller than 8 A. In contrast, the percentage increases to 80% and 81%, respectively, in membranes above the upper bound. In addition, only 29% of membranes below the upper bound possess a VF smaller than 0.5, in contrast to the 61% of structures above the upper bound. Finally, 32% of membranes below the upper bound have a density between 1 and 1.5 g cm À3 , in contrast to 64% of membranes above the upper bound. These results suggest that the structureperformance relationships described in this study are signicant and may provide target structural properties for designing future membranes with excellent separation performance.
Calculations of the heats of adsorption (Q 0 st ) can provide useful insight into a MOF's affinity for different gas molecules. These can be directly linked to the strength of host-guest interactions which are related to structural properties such as pore diameter (PLD, LCD) and porosity (VF). The Q 0 st distributions for the MOFs above and below the 2008 Robeson upper bound are presented in Fig. 8. For CH 4 , CO 2 and H 2 S, the distributions representing MOFs above the upper bound (red) lie further to the le than those below the upper bound (blue), implying that, in general, there are stronger host-guest interactions present in the MOFs above the upper bound. This is also reected in the average values of K 0 CH4 , K 0 CO2 and K 0 H2S which were a factor of 3.03, 3.90 and 4.73 larger for MOFs above the upper bound, respectively. From these results, we expect to see a reduction in the rate of diffusion for each biogas component. For the MOFs below the upper bound, we found that the average D 0 CH4 , D 0 CO2 and D 0 H2S were a factor of 1.87, 1.42 and 1.64 times larger than the MOFs above the upper bound. As permeability is calculated from the product of the solubility and diffusion properties of the material, we can see that the changes to K 0 i outweigh the changes to D 0 i and therefore, permeability and selectivity become larger for the MOFs above the upper bound.

Conclusions
In this work, a multi-stage screening methodology was designed to identify the most promising porous MOFs for biogas upgrading. The Cambridge structural database was pre-screened using a set of geometric criteria in order to remove structures with non-accessible surface area and pore diameters smaller than the individual components of biogas. The initial phase of screening predicted the single-component permeation of CO 2 , CH 4 and H 2 S from adsorption and diffusion calculations at innite dilution. In addition, Henry constants for the adsorption of H 2 O were calculated for each structure to determine the relative levels of hydrophobicity. Only the top hydrophobic MOF structures which exceed the Robeson upper bound (these are assumed to be the least affected by coadsorption of H 2 O molecules within the biogas stream) were considered for the nal phase of screening. For these materials, mixture selectivity and permeability at 298 K and 10 bar were predicted from equimolar CO 2 /CH 4 and H 2 S/CH 4 simulations. The top 8 MOFs identied in this work possess selectivity and permeability that rivals that of polymer-, zeolite-, ZIF-and MMM-based membranes. For the top candidates, the predicted values of CO 2 permeability range between 3.71 Â 10 4 barrer and 1.09 Â 10 5 barrer, CO 2 selectivity between 4.93 and 11.9, H 2 S permeability between 5.69 Â 10 3 barrer and 1.03 Â 10 5 barrer and H 2 S selectivity between 32.7 and 123 at the operating conditions of 10 bar and 298 K. In addition, the structural properties of the 1183 hydrophobic membranes with membrane selectivity values above the Robeson limit were examined. Structural property histograms indicated that MOF membranes with PLDs between 3.8 A and 6 A, LCDs between 4 A and 8 A, ASAs less than 1250 m 2 g À1 , VF between 0.3 and 0.5 and density between 1 g cm À3 and 1.5 g cm À3 appear the most frequently.
To verify the signicance of these properties, a direct comparison was made against the structural property distributions from low-performing materials found below the 2008 Robeson bound. It was found that 42% of low-performing materials possess PLDs and LCDs smaller than 6 A, and 8 A, respectively, which increases to 80% and 81% for membranes above the upper bound. In addition, 29% of low-performing materials possess a VF smaller than 0.5, in contrast to the 61% of membranes above the upper bound. A similar result was also noted when comparing the density distributions. This suggests that the structural properties that are frequently observed in the membranes above the upper bound are signicant and contribute to the superior separation performance of these membranes. Therefore, these results may inuence the future design of novel MOFs for biogas upgrading.

Conflicts of interest
There are no conicts to declare.