As a major component of greenhouse gases, carbon dioxide is a main cause of global warming [1-6]. Zeolitic imidazolate frameworks are promising porous materials for CO₂ capture because of their large accessible volumes, tunable functional groups, and relatively high thermal and chemical stability [7-14]. In particular, the properties of ZIFs are essentially determined by their network topologies [15, 16]. To this end, many computational efforts have been made to build a number of hypothetical ZIF models from specific topologies, from which candidates with the most promising CO₂ capture performance may be identified [17-20]. For instance, Smit et al. studied 1000 hypothetical ZIFs for CO₂ capture and identified 50 of them with the smallest parasitic energy as candidate materials [21]; Calero et al. constructed 197 hypothetical ZIF models from known zeolites and identified 10 candidates with the best Henry's constant ratio or adsorption selectivity [19]. Le et al. built 18 hypothetical ZIF models from 3 network topologies and identified two of them as potential candidates for CO₂ capture [20]. Despite all the success, the number of explored ZIF models is rather limited in comparison with other MOF materials [22, 23].

Here, we report a high-throughput study in model-building and screening of new ZIF structures for CO₂ capture from flue gas. We built 169, 898 ZIF models from 84, 949 hypothetical zeolite networks. These ZIF models were further evaluated by lattice energies, accessible volume to CO₂, the isosteric adsorption heat (Q_st) of H₂O, the ratio of Henry's constant (S_KH) of CO₂/N₂, percent regenerability (R%), CO₂ working capacity (ΔN_CO2), CO₂/N₂ adsorption selectivity (S_CO2/N₂) and adsorbent performance score (APS). Among such a large number of ZIF models we explored, we were able to identify 49 hydrophobic ZIF structures that might be superior to experimentally realised ZIF structures built from the same imidazolate linkers for CO₂ capture (Scheme 1).

The hypothetical zeolite networks we began with were obtained in our previous study, including 84, 292 ABC-6 networks generated by structure enumeration [24] and 657 tetrahedral networks generated by automated atom assembly [25, 26]. We wrote a perl script to build 169, 898 ZIF models automatically from zeolite networks by replacing the vertex Si atoms and the bridging O atoms in zeolite networks by Zn or Cu cations and imidazolate linkers, respectively. The ZIF models derived from ABC-6 networks are named ABC-ZIFs hereafter and those from automated atom assembly are named AAA-ZIFs. Then, the unit cell dimensions of the original zeolite networks were expanded by 1.9 times to approximate those of ZIF structures. All ZIF models were further optimised using the UNIVERSAL force field [27]. The highest allowed symmetries were imposed to these models throughout geometry optimisation. More details are shown in Table S1 (Supporting information). On the other hand, we optimised 24 Zn-ZIF and 5 Cu-ZIF structures that were experimentally realised [22]. All these calculations were conducted using Materials Studio software package [28]. Considering the fact that the hypothetical ZIF models with high lattice energies would be unstable and difficult to realise in experiment, we set the upper limits for the lattice energies of hypothetical ZIFs as the highest lattice energies observed in realised ZIFs. By doing this, we narrowed down the list of stable ZIFs to 83, 398 (Fig. S1 in Supporting information).

We calculated the pore-limiting diameters (PLDs) and the largest-cavity diameters (LCDs) for all stable hypothetical ZIFs using program Zeo++ [29]. The results are shown in Fig. S2a (Supporting information). Because all our ZIF models were built from zeolite networks, most of them possess large void spaces. In comparison to AAA-ZIFs (PLD: 2.5-7.5 Å; LCD: 5-10 Å), ABC-ZIFs exhibit larger void space for guest species (PLD: 4-7 Å; LCD: 10.5-14 Å) because of their featured cage structures. Then, the CO₂ accessible volumes and surface areas for these ZIFs were calculated using a spherical probe with a diameter of 3.3 Å, which corresponded to the kinetic diameter of CO₂ (Fig. S2b in Supporting information). 323 ZIF models were removed because of the absence of pores accessible to CO₂. We also calculated the CO₂ accessible volumes for experimentally realised ZIFs. Half of the realised Zn-ZIFs (12/24) and most of the realised Cu-ZIFs (4/5) were removed because of the absence of accessible pores.

To evaluate the CO₂ capture performance of these ZIF models, we conducted single-component Grand Canonical Monte Carlo (GCMC) simulations for the remaining 83, 075 ZIF models at 298 K. All simulations were performed using Materials Studio software package. All force field parameters were taken from literature [30-33]. Detailed information can be found in Table S2 (Supporting information). First, we calculated the isosteric adsorption heat, Q_st, and the Henry's constant, K_{H, H2O} of H₂O molecules in these ZIF models at infinite dilution. The results are shown in Fig. S3 (Supporting information). Only five of the remaining hypothetical ZIF models exhibited isosteric adsorption heat greater than 45 kJ/mol, and over 80% of them less than 30 kJ/mol. As expected, most of our ZIF models were not hydrophilic, indicating their potential applications under real flue gas environment with certain humidity [34].

Then, we separately calculated the adsorption of CO₂ and N₂ molecules of the remaining 83, 070 ZIF models and obtained their Henry's constant, i.e., KH, CO₂ and K_{H, N2} respectively. The ratio of Henry's constant of CO₂/N₂ was calculated using the formula:

(1)

which represented the ideal selectivity of the adsorbent at infinite dilution. The ratio of Henry's constant would be a good approximation of the true separation capability of the adsorbent under practical conditions, and its computational cost is much lower than regular multi-component GCMC simulations [35, 36]. According to our calculations, the Henry's constant ratios of the 83, 070 ZIF structures ranged from 5.83 to 8525.37 (Fig. S4 in Supporting information). We kept the top 10% models with the highest S_K for further simulations. Note that two of the realised ZIF structures were removed at this step because of the low Henry's constant ratios.

Next, we performed multi-component GCMC simulations for the remaining 8313 hypothetical ZIF models and 11 already-realised structures with the highest S_KH in CO₂/N₂/H₂O ternary systems. The simulations were conducted at 10 bar and 1 bar, corresponding to the gas adsorption pressure and the desorption pressure, respectively. The CO₂:N₂ molar ratio was set to 1:9, and a relative humidity of 80% was maintained in correspondence with the content of flue gas [37, 38]. Fig. 1 shows the density distribution of CO₂ molecules in two selected ZIF models during GCMC simulations in CO₂/N₂/H₂O ternary systems. The CO₂ molecules are mainly located close to the imidazolate linkers. Such affinity can also be observed in the molecular dynamics trajectories of individual CO₂ molecules (Fig. S5 in Supporting information).

Using the simulation date obtained from CO₂/N₂/H₂O ternary systems, we calculated the percent regenerability (R%), which is defined as the ratio of the working capacity (ΔN_CO2) and the adsorption capacity (N_CO2^ads) at high pressure:

(2)

(3)

where N_CO2^ads and N_CO2^des are loadings at adsorption and desorption conditions, respectively. R% represents the fraction of the adsorption sites that are regenerable at low pressure [39, 40]. The results are shown in Fig. S6 (Supporting information). 7959 of our remaining 8313 ZIF models exhibited R% ≥ 80%. In comparison, only 2 of the already-realised ZIF structures exhibited R% ≥ 80%. Because materials with low R% would not be suitable for practical cyclic utilization [39-41], we screened out 354 hypothetical models with R% < 80% and kept 7959 models for further evaluation. All the remaining high-R% ZIFs exhibited CO₂ isosteric adsorption heat less than 35 kJ/mol (Fig. S6), so the adsorbed CO₂ molecules could be removed easily at low pressure. This might be due to the highly porous nature of our ZIF models that were built from open-framework zeolite networks.

To identify the best structures for CO₂ capture from the remaining ZIFs, we further calculated adsorption selectivity (S_CO2/N2) and APS as follows [37, 39-42]:

(4)

(5)

where N_CO2^ads and N_N2^des are loadings at adsorption and desorption conditions, respectively; y is the mole fraction in the flue gas phase. All these data were calculated from GCMC simulations in CO₂/N₂/H₂O ternary systems.

The results are shown in Fig. 2. The maximum adsorption working capacity ΔN_CO2 of our hypothetical ZIF structures reached 2.37 mmol/g, which was much higher than that for experimentally realised structures (1.61 mmol/g). The highest adsorption selectivity S_CO2/N2 for our hypothetical ZIF models was 31.83, and 17.15 for already-realised ZIFs. The best APS for hypothetical and experimentally realised ZIFs were 68.70 mmol/g and 27.03 mmol/g, respectively. In general, the performance of the remaining hypothetical ZIFs was significantly better than already-realised ones. By ranking their adsorption working capacity, adsorption selectivity, and adsorbent performance score, we identified 49 hypothetical ZIF models that were among the top 100 in all categories. Detailed information of the 49 structures are listed in Table S3 (Supporting information). The isosteric CO₂ adsorption heat for the 49 hypothetical ZIF models were within the range of 24.26-30.99 kJ/mol, which were close to the ideal values reported by Snurr et al. [37].

In general, AAA-ZIF models exhibited better CO₂ capture performance than ABC-ZIFs, especially in adsorption selectivity and adsorbent performance score. This might be because the AAA-ZIF models possessed smaller void space than ABC-ZIFs, and relatively smaller pore openings are usually crucial for selective adsorption [4, 37]. Although our hypothetical ZIF models did not exhibit the best CO₂ selectivity among all types of ZIF materials (Table S4 in Supporting information), the results we obtained here were quite acceptable because we did not include any additional functional groups in our ZIF models, which could dramatically increase the affinity to CO₂ molecules [4, 9, 16, 43, 44].

Because of the lack of open metal site and amino group, the interactions between CO₂ molecules and our ZIF models were mainly physical adsorptions. To further understand the contributions of van der Waals interactions and electrostatic interactions in adsorptions, we performed single-component GCMC simulations separately with and without charge assignment, and calculated the corresponding adsorption isotherms for two selected ZIF models with the best working capacity and selectivity at 298 K [44]. The results are shown in Fig. S7 (Supporting information). For N₂ adsorption, electrostatic interactions almost had no contribution to the loading and isosteric heat for both ZIFs. For CO₂ adsorption, van der Waals interactions still dominated. However, electrostatic interactions did influence CO₂ adsorption capacity and isosteric heat, because CO₂ possessed a higher quadrupolar moment than N₂ [45, 46].

In conclusion, through high-throughput calculations, we constructed 169, 898 hypothetical ZIF models and identified the most promising candidates for CO₂ capture from flue gas by step-wise screening. Our hypothetical ZIF models exhibited much better performance than already-realised ZIFs built from the same imidazolate linkers. Our approach enables efficient construction and evaluation of such a large number of ZIF models that have never been achieved before. This will accelerate the experimental discovery of novel ZIF materials with properties superior to existing ones. In this work, we have only considered imidazolate linkers. In future study, we will add functional groups with different affinities to CO₂ to further tune the CO₂ capture performance.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (No. 2016YFB0701100), the National Natural Science Foundation of China (Nos. 21622102 and 21621001), the National 111 Project (No. B17020), Program for JLUSTIRT, and High Performance Computing Center of Jilin University.

Appendix A. Supplementary data

Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.cclet.2019.04.025.