Transferability of CO2 Force Fields for the Prediction of Adsorption

structures in an entire database. However, all .... sphere diameter (Di) retrieved from the IZA database,17 which suggest that zeolites with more conf...
4 downloads 5 Views 959KB Size
Subscriber access provided by University of Sussex Library

C: Surfaces, Interfaces, Porous Materials, and Catalysis 2

Transferability of CO Force Fields for the Prediction of Adsorption Properties in All-Silica Zeolites Jian Ren Lim, Chi-Ta Yang, Jihan Kim, and Li-Chiang Lin J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b02208 • Publication Date (Web): 19 Apr 2018 Downloaded from http://pubs.acs.org on April 19, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Transferability of CO2 Force Fields for the Prediction of Adsorption Properties in All-Silica Zeolites Jian Ren Lim,a Chi-Ta Yang, a Jihan Kim,b Li-Chiang Lina*

a

William G. Lowrie Department of Chemical and Biomolecular Engineering, The Ohio

State University; 151 W. Woodruff Ave., Columbus, Ohio 43210, USA b

Department of Chemical and Biomolecular Engineering, Korea Advanced Institute of

Science and Technology; 291 Daehak-ro, Yuseong-gu, Daejeon 305-338, Republic of Korea

* Corresponding author: L.-C.L ([email protected])

Abstract We present a systematic and comprehensive investigation of available CO2 force fields for their predictions in adsorption properties in 156 geometrically diverse zeolite structures. The comparison reveals that a large discrepancy in the predicted properties by more than two orders of magnitude may exist. Especially, variation predicted by different force fields appear to be more pronounced for zeolites with more confined pore features, which can be attributed to the repulsive characteristics of force fields. The discrepancy especially impacts zeolites that are deemed to be the best materials for CCS, indicating that the predictions on the best materials can drastically differ based on the choice of force fields. To develop accurate and fully transferable force fields, in this work, we show that the inclusion of adsorption uptake at a highpressure region (or saturation loading) as well as the diffusion coefficient can be of utmost importance. These properties can be used as indicators to the repulsive behaviors between gas molecules and the framework. Mixture isotherms have also been identified to be potentially useful for the same purpose. Moreover, we have also demonstrated that interaction energies computed using ab initio methods can be useful references to ensure a newly developed force field is capable of describing the energy surface at an atomic level. Overall, the outcomes of this study will be instrumental to the future development of accurate and transferrable force fields, which is critical for future large-scale computational studies.

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 30

1. Introduction Since the late twentieth century, anthropogenic carbon dioxide pollution (i.e., global warming and climate change) has become one of the most challenging issues.1,2 Due to the continuing use of fossil fuels, such challenge is projected to become even more pronounced in the future. To reduce the emissions of CO2 into the atmosphere, carbon capture and sequestration (CCS)3–5 would represent a viable short-term solution except for its high-energy costs. Amine scrubbing method is currently the state-of-the-art approach to capture CO2.6 However, due to the chemical bonding between CO2 and amine molecules, a large amount of energy is required for the regeneration of amine sorbents. As such, developing more energyefficient approaches to capture carbon dioxide is of critical importance. To this end, nanoporous materials have been regarded as promising candidates.7–12 These materials possess a variety of desirable properties (e.g. large surface area, moderate heat of adsorption, low heat capacity) to efficiently separate CO2 molecules from the flue gas mixture. Specifically, siliceous zeolites have drawn particular attention due to their thermal stability, chemical robustness, and their hydrophobic nature.13,14

Microporous zeolites are highly tunable with a wide range of diverse geometrical features. Due to their tunability, more than 100,000 hypothetical zeolite structures have been theoretically predicted,15,16 although only approximately 230 zeolites have been synthesized as included in the International Zeolite Association database17 (denote as IZA zeolites hereafter). Such a large material space provides great opportunities in discovering optimal adsorbents. At the same time, it also makes the materials search unprecedentedly difficult. To identify promising adsorbents, computational approach can play an important role. By using state-ofthe-art molecular simulations with the modern computer architecture, a large number of structures can be investigated in an efficient and effective manner. As an example, recently, a number of large-scale computational screening studies have been deployed to explore over one hundred thousand zeolite materials for carbon capture applications.12,18–22 One of our recent studies has identified optimal zeolite structures that may cut the energy-requirement for CCS by as much as 35% as compared to that of the amine method.12 Similarly, zeolite materials have also been computationally explored as adsorbents in other applications such as paraffin/olefin separations,23 natural gas purifications,24–29,30 Kr/Xe separations.31–33.

ACS Paragon Plus Environment

2

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A key component in these large-scale computational studies is the accuracy of the adopted force fields to describe intermolecular interactions, thus providing sensible predictions of the gaseous adsorption properties in the structures of interest. To the best of our knowledge, there are at least six available force fields34–39 reported in the literature that have been shown to provide good predictions for CO2 adsorption in siliceous zeolites. The one developed by GarcíaPérez et al.36 appears to be most widely used, possibly due to its early development (i.e., developed 10 years ago). For simplicity and practicality while carrying our large-scale studies, a critical assumption made in all of the screening studies to date was that the adopted potential can be fully transferable; a particular set of parameters can be applied to calculate properties for structures in an entire database. However, all of these available force fields have been parameterized and/or validated based with only a few structures. As a result, it remains unclear whether such a transferable assumption would hold for structures that are diversified in structural features. Additionally, it remains to be seen if the predictions of CO2 properties as well as the evaluation of zeolites’ performances in CCS would differ drastically when different sets of force fields are used. Accordingly, in this study, we have carried out a systematic computational study to compare the predictions made by available force fields for CO2 adsorbed in zeolites with diverse structural features. A large-scale computation was conducted to calculate CO2 adsorption properties in 156 all-silica zeolites using each of the available force fields. Insights toward the variation in predictions made by different force fields were achieved. The observed discrepancy in computed properties was also correlated to the structural features and the energy caracteristics of force fields. Furthermore, several key properties such as saturation loading and diffusion coefficient were proposed to be included for the validation of these force fields, and more importantly, for the future development of a fully transferable force field. Additionally, dispersion-corrected density functional theory calculations were employed to compute the interaction energies of CO2 molecules with zeolite frameworks and compared with the prediction of the available force fields. This information was also identified to be critically important towards the development of transferable potentials.

2. Computational details Monte Carlo simulations, implemented in the open-source RASPA 2.0 package,40,41 were carried out in this study to compute gaseous adsorption properties of CO2 in all-silica

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 30

zeolites materials. A total of 156 structures from the IZA database17 with diverse structural features were investigated (i.e., structures accessible to CO2 with an orthorhombic or hexagonal unit cells). As noted above, six available force fields, including that developed by Vujic et al.,34 Bai et al.,35 García-Pérez et al.,36 Cygan et al. (so-called CLAYFF),37 García-Sánchez et al.,38 and Fang et al. (so-called D2FF)39 were used in this study to predict CO2 adsorption properties in zeolites. To describe intermolecular interactions, both van der Waals (6-12 Lennard-Jones) and long-range Coulombic potentials were adopted in all of these force fields. The parameters of these force fields are summarized in Table S1 of the Supporting Information (SI). All of these force fields except that of Vujic et al.34 and CLAYFF37 were developed based on rigid models (i.e., assumed rigid zeolite frameworks). Interestingly, the results of Vujic et al. indicated that their proposed potential might also be used with rigid assumption given that both the flexible and rigid models resulted in nearly no difference in the calculated adsorption properties. Meanwhile, CLAYFF was in fact initially not developed for CO2 adsorption in clay materials but to simulate adsorption of water. However, as shown by Fang et al.,39 CLAYFF (i.e., in combination with the EPM242 force field for CO2 molecules with a rigid assumption) may also be used to simulate CO2 adsorption in zeolites, and the prediction is in a decent agreement with that measured experimentally. As a result, we included both of the aforementioned force fields in this study and compared all of these six models with the assumption of rigid frameworks. In order to evaluate the carbon capture performance of these studied zeolites under a simplified flue gas condition (i.e., CO2 and N2 mixture), we also computed N2 adsorption properties in these structures. For this purpose, the N2 force field developed by García-Pérez et al.36 was used.

A variety of adsorption properties were computed in this work, including pure component and mixture adsorption isotherms, Henry coefficients, and the free energy landscape of CO2 adsorption inside the zeolite framework. Monte Carlo simulations in a grand canonical ensemble (GCMC)43 were used to compute the adsorption isotherms (for both pure component and mixture), while the Widom particle insertion method44 was employed to calculate the Henry coefficient of adsorption. Millions of configurations were sampled via a variety of Monte Carlo moves to obtain statistically sound averages. Monte Carlo simulations in a canonical ensemble (NVT) were also conducted to probe the free energy landscape of CO2 adsorbed in the framework. The carbon positions of CO2 molecules were analyzed and projected onto a twodimensional (2D) plane along a certain axis. For better presentation, the free energy landscape projected onto a 2D plane was presented in the form of density maps. Warmer colors indicate a

ACS Paragon Plus Environment

4

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

higher probability of finding CO2, indicating a lower free energy. In all of these calculations, the 6-12 L-J potential was truncated and shifted to a cutoff radius of 12 Å, while the long range Coulombic interactions were computed by the Ewald summation technique.41,45 The simulation box of a given studied structure was composed of one to multiple unit cells, ensuring that its dimension is at least twice the cut-off radius along each direction for periodicity. We note that, in these calculations, blocking of structure was not considered. Additionally, density functional theory calculations implemented in Quantum Espresso46 code with a planewave basis were also performed to compute interaction energies between a single CO2 molecule and a zeolite framework for a number of randomly inserted configurations. The non-local vdW-DF2 exchangecorrelation functional47 was adopted for dispersion energies. In these calculations, TroullierMartins type norm-conserving pseudopotentials48 were utilized, while kinetic energy and charge density cutoffs were set to be 80 Rydberg and 320 Rydberg, respectively.

3. Results and discussions In this section, we first show the comparison between available force fields for their predicted CO2 adsorption properties. A considerable discrepancy as high as several orders of magnitude was found, and such discrepancy was also identified to significantly influence the evaluation of zeolites’ performance in CCS. We have demonstrated a strong correlation between the discrepancy and structural properties of zeolites such as the largest included sphere diameter (Di) retrieved from the IZA database,17 which suggest that zeolites with more confining features may potentially lead to large sensitivity in the predictions. To facilitate the validation of force fields and the future development of fully transferable force fields, several properties such as saturation loading and diffusion coefficient are proposed as key references. We have also carried out density functional theory calculations to compute the interactions between CO2 and zeolite frameworks, serving as references to potentially probe the accuracy of these available force fields.The references for models of Vujic et al., Bai et al., García-Pérez et al., CLAYFF, García-Sánchez et al., and D2FF discussed in the following sections can be found in Ref [34 – 39], respectively.

3.1. Predictions by different force fields and their impacts Henry coefficient,44 the initial slope of adsorption isotherms, represents the adsorption strength of CO2 in the low-pressure region where the intermolecular interactions between host

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 30

(zeolites) and guest (CO2) are dominant. Given all of these available force fields were developed to only parameterize the host-guest interactions with guest-guest interaction parameters adopted directly from widely used CO2 molecular models such as TraPPE49 and EPM242, the Henry coefficient of adsorption can serve as an important probe and the first step to compare these available force fields. It should also be noted that Henry coefficients have been found to highly correlated with the separation performance of adsorbent materials in CCS.11,12 Figure 1 summarized the Henry coefficients of CO2 in all studied zeolites predicted by different force fields. In this figure, for clarity, we have chosen to present the data in terms of the ratio of the Henry coefficients with respect to the values predicted by D2FF. The order of zeolites were sorted based on the ratio predicted by the model of García-Pérez et al. Also, to serve as the guide to the eyes, two dashed lines were given to indicate ratios of 0.5 and 1.5 (corresponding to a deviation of -50% and +50%, respectively), and we propose that predictions within such a range should be deemed as acceptable. As is well known, the atomic structures of a given framework type of zeolites determined from experiments may be slightly different, and such a difference may lead to variations in computed predictions.50 Using MFI, one of the most well-studied zeolite structures, as an example, we carried out molecular simulations using different available atomic structures of MFI to compute their adsorption properties. As shown in SI Figure S1, our results indeed show that there may be a variation of as high as 50% due to the structure difference. With such an acceptable range as the guide, Figure 1 shows that the predictions made by different force fields evidently do not agree with each other and the differences can be surprisingly more than two orders of magnitudes. The model of GarcíaSánchez et al. appear to be distinct from the others, as its predictions for most sturctures are out of the range of +/- 50% w.r.t. to that by D2FF. Although the predictions made by the other five force fields seem to agree reasonably well (i.e., within +/- 50%), a non-negligible portion (i.e., approximately 40%) of zeolites studied herein still show a notable prediction variation .

ACS Paragon Plus Environment

6

Page 7 of 30

10.0

Henry coefficient ratio (w.r.t. D2FF)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

5.0

MFI CHA

1.5 1.0

0.5 Vujic et al. Bai et al. García-Pérez et al. CLAYFF García-Sánchez et al. 0.1

IZA structures Figure 1. Ratio of CO2 Henry coefficients predicted by various force fields w.r.t D2FF in 156 zeolites at 300 K. The order of zeolite structures is sorted based on the ratio predicted by García-Pérez et al.36 Two dashed lines indicate a deviation of +/-50% from unity to serve as the guide to the eyes. Two most widely studied zeolites, MFI and CHA, are highlighted in the figure. We note that this figure is focused in the region between a henry coefficient ratio of 0.05 and 10 to intensify the comparison between force fields (see SI Figure S2 for a full-scale plot). The same plot in terms of absolute values instead can be seen in SI Figure S3. Additionally, similar ratio plot w.r.t. García-Pérez et al. with the order of zeolites sorted by D2FF can be found in SI Figure S4.

It appears that one may categorize force fields into several groups. Besides the distinct characteristic of the force field model by García-Sánchez et al., the results shown in Figure 1 suggest that the predictions made by models of García-Pérez et al. and Vujic et al. lead to similar results, while CLAYFF and D2FF are alike to each other. Results of the model by Bai et al. also resemble the predictions of CLAYFF and D2FF well but appear to consistently estimate a weaker adsorption strength. To quantify such an observation, we have also computed the Spearman rank correlation51,52 between different force fields. The obtained correlation rank between two sets of data represents their similarity from the range of -1.00 (inversely

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 30

proportional), 0.00 (no correlation) to 1.00 (directly proportional). Accordingly, the Spearman rank correlation computed using the eq. (1) shown below indicates a one-to-one relationship between force fields. 6 ∑ 

 = 1 − 1 − 1 Where  is the square of difference between two ranks, and n is the number of samples. The Spearman ranks between available force fields are summarized in Table 1. It is observed that the ranks can vary from 0.17 (essentially not correlated) to 0.95 (highly correlated). The model by García-Sánchez et al. is fundamentally different from others (i.e., with a coefficient of as low as 0.17 w.r.t. to other force fields). The similarity between models of García-Pérez et al. and Vuijc et al. can be found from the high correlation value of 0.91. Meanwhile, strong correlations can also be seen between D2FF, CLAYFF, and Bai et al (0.88-0.95). Although the model by Bai et al. consistently predicts a lower Henry coefficient compare to that of D2FF and CLAYFF, the relative strength amongst studied materials remain roughly the same. Overall, these results quantitatively confirm the observations made from Figure 1. Interestingly, in contrast to the adsorption properties at the infinite diluted condition, the Spearman rank coefficients for uptakes predicted at a high pressure (i.e.1 GPa ) shows a diminished discrepancy as can be seen from SI Table S2. This appears to be reasonable as the uptake at a higher pressure region will also be largely dominated by the structural features, making the difference to be relatively smaller.

Table 1. Spearman rank correlation of Henry coefficients at 300 K between different force fields. Vujic et 34

Vujic et al.

34

35

Bai et al.

García-Pérez et al. CLAYFF

36

37

Bai et 35

GarcíaPérez et 36 al.

CLAYFF

37

GarcíaSánchez 38 et al.

D2FF

39

al.

al.

1.00

0.74

0.92

0.90

0.30

0.81

0.74

1.00

0.65

0.87

0.73

0.89

0.92

0.65

1.00

0.88

0.15

0.80

0.90

0.87

0.88

1.00

0.38

0.95

0.30

0.73

0.15

0.38

1.00

0.48

0.81

0.89

0.80

0.95

0.48

1.00

38

García-Sánchez et al. D2FF

39

ACS Paragon Plus Environment

8

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The development of these force fields generally involves the fitting of force field parameters to experimental adsorption isotherms in a set of structures (so-called training set), and only few of them were parameterized based on reference interactions derived from quantum chemical calculations (e.g., D2FF). Another set of structures (denoted as the validation set) will be generally subsequently used to validate the parameterized force fields, thus rendering a transferable force field. It should be noted that the forcie field of García-Sánchez et al. was parameteried to model CO2 adsorption in aluminosilicate zeolites instead of all-silica ones, but the authors have also demonstrated that the developed potential can be applicable as well to siliceous zeolites. While most of these force fields have claimed their transferability, our results indicate that a significant discrepancy does exist. Table 2 summarizes the adsorption properties of zeolites that were included in the development of the force fields studied herein. Interestingly, the predictions made by different force fields in these structures were found to yield similar results. For example, MFI and CHA, two of the most widely studied zeolites, were essentially included in the parameterization and/or validation of all of these available force fields due to their wide availability of experimental data. As shown in Table 2 and Figure 1, the adsorption properties in MFI and CHA are insensitive to the choice of force fields. This indicates that the large discrepancy observed in the predictions made by various potential may be likely attributed to the lack of diversified selections of structures in the force field development. Most of the included structures generally have a fairly wide cages and/or channel with a diameter of ~ 5.5 Å or larger. This also implies these available potentials have not yet been validated against any experiment data for zeolites which possess a relatively smaller confinement. A detailed discussion can be seen below, which indicates that structures with more confining pore features indeed lead to a larger discrepancy.

Table 2. Ratio of Henry coefficients w.r.t D2FF at 300 K for zeolites that have been used in the validation and parameterization of the six studied force fields. Zeolites that were used in parameterization and validation sets were labelled with superscript ‘a’ and ‘b’ respectively. The largest included sphere diameters (Di, in Å) of structures are listed in parentheses. The absolute values of Henry coefficients for these materials can be found in SI Table S3 in SI) Henry Coefficient Ratio (w.r.t. D2FF) Vujic et

Bai et

García-

al.34

al.35

Pérez et al.36

GarcíaCLAYFF37

ACS Paragon Plus Environment

Sánchez et al.38

D2FF39

9

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 30

MFI (6.36)

0.92ab

0.51ab

0.79ab

0.93

1.21b

1.00b

FER (6.31)

1.23ab

0.67

0.73

0.98

3.23

1.00

CHA (7.37)

ab

0.95

0.57

b

0.85

0.97

1.44

1.00a

BEA (6.68)

0.92ab

0.76

0.64

1.01

2.26

1.00

LTA (11.05)

1.21a

0.65

0.66

0.95

3.31

1.00

DDR (7.66)

1.95

0.62b

1.16

1.04

2.19

1.00b

MOR (6.70)

1.23

0.26

2.29b

1.19

0.31

1.00

As the goal of computational screenings is to identify promising candidate materials, at this point, it is of utmost importance to unravel the influence of the discrepancy in adsorption predictions on the evaluation (ranking) of materials for CCS. To evaluate zeolites as adsorbents to capture carbon,

a metric denoted as parasitic energy11,12 was previously proposed to

compute the energy penalty per amount of captured CO2 imposed on power-plants when an adsorbent material is used. This metric provides a comprehensive evaluation of materials, which takes the adsorption, desorption, and compression steps involved in CCS into account. Based on the properties predicted using different force fields, the parasitic energy values for each zeolite is summarized in Figure 2. Our results indicate that the identified top-ranked materials may drastically differ based on the choice of force fields due to the discrepancy in the prediction of CO2 adsorption properties. Specifically, a zeolite predicted as an optimal material (i.e. lower parasitic energy) using one set of a force field might not be as promising according to another force field. Table 3 shows the ten top-ranked structures evaluated by the model of García-Pérez et al. and the corresponding rankings by other five force fields. Although several of these ten top-ranked candidates are still amongst the top-ranked materials using other force fields, the dissimilarity can be obviously acknowledged. For example, the first-ranked candidate material, AWW, predicted by the model of García-Pérez et al. is predicted to be 47th by D2FF.

ACS Paragon Plus Environment

10

Page 11 of 30

8x103 Vujic et al. Bai et al. García-Pérez et al.

Parasitic energy (kJ/kg CO2)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

CLAYFF García-Sánchez et al. D2FF

6x103

4x103

2x103

0

IZA structures

Figure 2. Parasitic energy predicted by different force fields. The order of zeolite structures is sorted based on the parasitic energies predicted by the model of García-Pérez et al. Table 3. Rank of zeolites based on parasitic energy. The largest included sphere diameter (Di) of each structure is listed in the parenthesis. The absolute parasitic energy of these ten structures shown in the table can be found in SI Table S4. Ranking of structures based on parasitic energy

AWW (7.48)

GarcíaPérez et al.36 1

Vujic et al.34

Bai et al.35

54

32

53

GarcíaSánchez et al.38 15

LOS (6.36)

2

9

20

28

7

56

EWT (7.75)

3

100

31

5

113

4

MAR (6.35)

4

18

16

45

4

43

WEI (4.19)

5

4

63

9

149

12

TOL (6.37)

6

12

27

32

13

46

JBW (4.29)

7

1

5

3

148

6

GIS (4.97)

8

2

2

2

96

1

WEN (5.53)

9

19

1

1

71

32

ATN (5.91)

10

3

3

6

25

5

37

CLAYFF

ACS Paragon Plus Environment

D2FF39 47

11

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 30

3.2. Underlying physics behind the observed discrepancy While a huge discrepancy was identified and discussed above, it remains of utmost importance to understand the relationship between the prediction difference and the structural features of zeolite. Although there is no indication at this point to suggest which force field is more reliable, such a relationship can point to a group of zeolites containing specific features that can show order-of-magnitude variations in predictions if different force fields are used. Extra care should then be given to these structures in screening studies. As is well known, the strength of CO2 adsorption

is largely correlated to the size of the confined space.7,53,54

Specifically, the largest included sphere diameter,55 Di, is the diameter of the largest sphere that can fit in a structure, serving as a representative index for representing the extent of structural confinement. Figure 3 shows the ratio of Henry coefficients (i.e., w.r.t. D2FF) as a function of Di. The results indicate that structures with smaller pore size appear to show significantly larger deviation in the predictions made by different force fields. When Di ≤5 Å, predictions in Henry coefficients can differ by more than two orders of magnitude. Within this region, both models by García-Pérez et al. and Vujic et al. tend to predict a notably stronger interaction compared to that of D2FF, whereas models of Bai et al. and García-Sánchez et al. (i.e., in particular GarcíaSánchez et al.) appear to yield a much lower adsorption strength. For Di > 5.5 Å, the predictions made by all force fields except that by the model of García-Sánchez et al. agree reasonably well with each other (mostly within +/- 50%), although the model of Bai et al. tends to systematically underestimate the Henry coefficients. Interestingly, while García-Sánchez et al. leads to a considerably smaller Henry coefficients when Di ≤ 5 Å, it predicts a relatively higher henry coefficient for structures with Di between 5 Å and 11 Å, again leading to a behavior distinct from five other models. Overall, we found that structures with a small pore (i.e., Di ≤ 5 Å) would require extra care in the choice of potential, where the predictions are more sensitive to force field parameters. However, as shown in Table 2, none of the structures used previously for the development and/or validation of force fields have a pore feature in this region. MFI and CHA have a Di of 6.36 Å and 7.37 Å, respectively, both of which are in the region that force field predictions are considerably less sensitive to the selection of force fields. The Henry coefficient ratios as a function of other structural features including density and accessible volume ratio are respectively shown in SI Figure S5 – S6. Similarly, these figures illustrate that structures that are more confined (e.g., larger density and smaller pore volume) tend to show a larger discrepancy in predictions.

ACS Paragon Plus Environment

12

Page 13 of 30

10.0

Henry coefficient ratio (w.r.t. D2FF)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Vujic et al. Bai et al. García-Pérez et al. CLAYFF García-Sánchez et al.

5.0

1.5 1.0

0.5

0.1 2

4

6

8

10

12

14

16

Di (Å)

Figure 3. The ratios of Henry coefficients at 300 K w.r.t. that of D2FF as a function of the largest included sphere diameter (Di). A larger prediction deviation is observed when Di is less than approximately 5 Å. Dotted lines represent ratios of 0.5 and 1.5. A full-scale plot can be seen in SI Figure S7.

The observed correlation discussed above is anticipated to closely relate to the characteristics of force fields, given the discrepancy in predictions should be attributed to difference in the host-guest potential energy curves of force fields. To probe this, we collected a set of ten thousand random CO2 configurations in zeolite DON, a structure possessing large one-dimensional (1D) channels, and calculated the energy and the shortest distance from the surface (i.e., C(CO2) and the framework atoms) of each of the configurations. A histogram analysis was then carried out with a distance bin size of 0.2 Å. For all configurations within the bin, we computed the Boltzmann-weighted average energy using the equation shown in eq. (2):

<  > =

∑    − /  2 ∑   − / 

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 30

Where U is the interaction energy,  is the Boltzmann constant, T is the absolute temperature, and N is the total number of configurations in the bin. The resulting potential curve of each force field is shown in Figure 4, which resembles the typical shape of Lennard-Jones potential. Surprisingly, the repulsive behavior of these potential curves shows vast differences, and the difference was found to be strongly correlated to the observed discrepancy in the prediction of Henry coefficients discussed above. For a distance of smaller than 3.0 Å, the model by GarcíaSánchez et al. predict a significantly stronger repulsive energy followed by that of Bai et al. Such a highly repulsive behavior in turn predicts a significantly lower Henry coefficients. In contrast, the model of García-Pérez et al. followed by that of Vujic et al. predict the least repulsive behavior, thus yielding larger Henry coefficients for structures with smaller pore sizes. At a larger distance (i.e. distance ≥ 3.2 Å), energies computed using all force fields except that of García-Sánchez et al. are similar, and the locations of the potential well are also in agreement with each other (i.e., at approximately 3.0 Å). The model of García-Sánchez et al. shows a drastically different behavior to have a notably deeper potential well at a larger distance of ~3.6 Å, which correspond to that distinct behavior in Henry coefficient predictions for structures with a Di between 5.5 and11 Å. Overall, a clear connection between the discrepancy in the predictions and the characteristic of the potential curves was established. The difference in the repulsive behavior of potentials was identified to be the key, thus making confined zeolites to be more sensitive to the choice of force fields.

ACS Paragon Plus Environment

14

Page 15 of 30

3x103

Boltzmann-Weighted Energy (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Vujic et al. Bai et al. García-Pérez et al. CLAYFF García-Sánchez et al. D2FF

2x103

1x103

0

-1x103

-2x103

2

3

4

5

6

Distance (Å)

Figure 4. Interaction energy between guest (CO2) and host (zeolite DON) predicted by various force fields as a function of the minimum distance between C(CO2) and the nearest atom of the host framework.

At this point, it is important to also investigate the role of Coulombic interactions in the observed discrepancy in adsorption predictions between different force fields. For this purpose, Henry coefficients with the electrostatic interactions turned off were computed and the results are shown in SI Figure S8. In comparison to Figure 3 (i.e., the data that includes both 6-12 Lennard-Jones and Coulombic interactions), the observed large discrepancy in the predicted adsorption properties remains similar. Orders of magnitude differences in the Henry coefficients are still observed, which are notable in particular for structures with more confining pore features (i.e., Di ≤ 5 Å). Therefore, this suggests that the prediction discrepancy between force fields is mainly driven by the van der Waals contributions, specifically the repulsive behaviors of different force fields as noted above. While the overall trend remains similar, we found that the Henry coefficients computed without Coulombic interactions by the models of García-Sánchez

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 30

et al., Vuijc et al., and Bai et al. show noticeable differences from the values calculated with

Coulombic

interactions.

For

instance,

without

including

the

Coulombic

contributions, the model of García-Sánchez et al. seems to predict an even larger discrepancy compared to other force fields for structures with Di > 5.5 Å. This result in fact appears reasonable as the atomic charges assigned by the model of GarcíaSánchez et al. are drastically distinct from the charges proposed in other force fields (see Table S1). Figure S9 further shows the comparison of the Coulombic contributions between force fields for each configuration in the abovementioned set of random CO2 positions in zeolite DON. As expected, the Coulombic interactions predicted by different force fields are generally similar and show a proportional relationship. However, owing to the difference in the assigned atomic charges, the model of García-Sánchez predicts the smallest Coulombic contribution in magnitude.

3.3. Development of transferable force fields The outcomes of the section 3.1 and 3.2 are instrumental to the potential validation, and more importantly, to the future development of transferable force fields. To determine the accurate force field models, in this study, we propose that several key references including, saturation loading (or uptakes at a high-pressure region), mixture component isotherm, and/or diffusion coefficient should be included. We will also demonstrate in this section that ab initioderived (e.g., DFT-derived) potential may result in a force field that not just reproduce ensemble averaged properties but also capture the detailed energy landscape of gaseous adsorption in nanoporous materials.

3.3.1 Saturation loading and mixture component isotherms One would expect that the force field models that predict relatively high repulsive would simultaneously predict relatively low saturation loadings due to smaller effective pore size. Similar to the analysis conducted for Henry coefficients discussed above, we have therefore also computed the saturation loading of CO2 in all of the zeolites studied herein using all force fields. In this analysis, for simplicity, the CO2 saturation loading in each structure is taken as the uptake value at a sufficiently large fugacity of 1 GPa to ensure saturation region is achieved. Figure 5 shows the ratio of the saturation loading w.r.t D2FF as a function of Di. In comparison

ACS Paragon Plus Environment

16

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

to D2FF, we found that models of García-Pérez et al. and Vujic et al. indeed predict a higher saturation loading, whereas the model of García-Sánchez et al. underestimates saturation loading the most. Consistent with Figure 3 (i.e., Henry coefficients vs. Di), a larger derivation in the predicted saturation value is observed for Di ≤ 5 Å. Interestingly, even for structures with a Di > 5.5 Å where Henry coefficient predictions by different force fields tend to be similar to all structures, a notable systematic deviation in the predicted saturation uptakes still exist. This suggests that a particular focus on the fitting to saturation loading (or high-pressure region values) in the development of force field may be critically important to ensure the repulsive characteristic of potential to be properly described. However, in the development of these available force fields, focuses have been primarily made on the fitting to adsorption uptakes in a relatively lower pressure region.

Figure 5. Saturation loading (i.e., uptake at a fugacity of 1 GPa at 300 K) predicted by each available force field as a function of largest included sphere diameter (Di).The results are presented by ratios with respect to that predicted by D2FF.

In addition, as one would expect from the ideal adsorption solution theory (IAST)56 for predicting mixture isotherms, the saturation loading (or high pressure region uptake) of each

ACS Paragon Plus Environment

17

The Journal of Physical Chemistry

individual component can influence its competitive nature57,58 with other gases. Accordingly, we hypothesize that the difference in the repulsive behavior of force fields can be also probed from the mixture component measurements, which can be therefore potentially be used as another reference in the force field development. As a proof of concept, zeolites with a consistent prediction in Henry coefficients by different force fields were studied for their CO2/N2 mixture isotherms. Figure 6 shows the mixture isotherms of CO2/N2 (i.e. 14% CO2 and 86% N2) in zeolite OBW, which clearly validates the aforementioned hypothesis. Due to the less repulsive nature of the model by García-Pérez et al. (i.e., a higher high-pressure loading), it predicts a more competitive CO2 adsorption over N2. As a result, as oppose to that of García-Sánchez et al. and D2FF, the N2 uptake is correspondingly predicted to decrease at a total fugacity of 50 bar or higher. Accordingly, having mixture isotherms may be potentially useful to verify if a physical behavior in repulsion has been properly captured. However, unlike the saturation loading or high-pressure region loading of pure component isotherms, mixture isotherms are may be too indirect to be used in the potential parameterization but rather an extra information for validation. Moreover, an additional uncertainty (in this case, uncertainty in the force field of N2) would also be involved. 8

Adsorption loading (mol/kg)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 30

6

García-Pérez et al. CO2 García-Sánchez et al. CO2 D2FF CO2 García-Pérez et al. N2 García-Sánchez et al. N2 D2FF N2

4

2

0 100

101

102

Fugacity (bar)

ACS Paragon Plus Environment

18

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 6. CO2/N2 (i.e. 14% CO2 and 86% N2) mixture isotherms in zeolite OBW predicted by using models of García-Pérez et al., García-Sánchez et al., and D2FF. As noted in the previous section, the adopted force field for N2 was taken from the work of García-Pérez et al.36

3.3.2. Diffusion coefficients As shown in Figure 2 and Table 2, the discrepancy in the Henry coefficients predicted by different force field generally becomes negligible for materials with less confined pore structures. For example, zeolite CHA has a relatively large Di of 7.37 Å and a reasonably consistent prediction in Henry coefficient (i.e., well within +/- 50%, see Table 2). However, although CHA has large cages, these cages are connected by a considerably smaller channel. The pore limited diameter of CHA, representing the diameter of the largest sphere that can transverse the structure, has a value of 3.72 Å. Therefore, the potential energies of CO2 molecules residing in such small channels may likely be very sensitive to the choice in the force field parameters. As a proof of concept, Figure 7 shows the density maps of carbon of CO2 projected onto the x-y plane in zeolite CHA computed by models of (a) D2FF, (b) García-Pérez et al., and (c) GarcíaSánchez et al., respectively. These density maps were obtained by analyzing millions of configurations from NVT simulations as introduced previously in computational details section. In the density maps, warmer colors indicate more favorable configuration of adsorption. These density maps demonstrate that the energy landscape of CO2 predicted by different force field can indeed vary substantially, although the resulting Henry coefficients are similar. The model of García-Pérez et al. predicts a very strong, favorable CO2-framework interaction inside the channels, whereas that of García-Sánchez et al. determines the channels to be nearly inaccessible to CO2. To further quantify the difference, Figure 7 (d-e) shows the corresponding Helmholtz free energies along the crystallographic “b” direction. Due to the discrepancy in interactions at atomic levels, the free energy barriers for transport predicted by different models can accordingly vary from 4.1 KBT (D2FF), 2.9 KBT (García-Pérez et al.), to 8.7 KBT (GarcíaSánchez et al.). Per the transition theory, more than two orders of magnitude difference in the CO2 diffusion properties predicted by these force fields can be expected. Therefore, experimental-determined diffusion coefficients can be an important reference in the future validation and parameterization of transferable force fields. It is important to also note that diffusion properties of gases in porous materials are critical in a variety of applications such as adsorption separations and membrane processes. Developing a force field that is capable of accurately predicting such properties is certainly of utmost importance.

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry

(a)

(c)

(b)

1.0 27.35

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.2

0

0

0

y (Å)

y (Å)

0.2

0

0

0

y (Å)

23.69

(f)

∆F (kBT) Coordinate (Å)

0.2

23.69

(e)

(d)

0.4

0

0

23.69

0.4

∆F (kBT)

0

x (Å)

1.0 27.35

x (Å)

x (Å)

27.35

∆F (kBT)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 30

Coordinate (Å)

Coordinate (Å)

Figure 7. (Upper panels, a-c) The density maps of C(CO2) adsorption projected onto the X-Y plane in zeolite CHA as predicted by (a) D2FF, (b) the model of García-Pérez et al., and (c) the model of García-Sánchez et al. respectively. (Lower panels, d-f) The respectively corresponding free energy diagram along the crystallographic “b” direction.

3.3.3 Reference energies from ab initio calculations To ensure that the energy predicted for each atomic configuration is reliable for validation and parameterization purposes, having interaction energies computed from ab initio quantum chemical methods such as density functional theory (DFT) calculations can also be important references to include in force field development. Using zeolite CHA again as an example, we have employed DFT calculations implemented in Quantum Espresso with nonlocal van der Waal (vdW-DF2) exchange-correlation functional to compute the interaction energies of CO2 with the zeolite framework for a set of 200 configurations obtained from an NVT simulation at 300 K. Figure 8 shows the comparison between DFT-calculated energies and force field-computed energies. Comparing to that predicted by models of García-Pérez et al. and García-Sánchez et al., the energies predicted by D2FF seemingly corresponds reasonably well with DFT calculations but with constant energy offset of approximately 1000 K. Considering

ACS Paragon Plus Environment

20

Page 21 of 30

the uncertainties involved in the choice of the exchange-correlation functional, the agreement can be deemed as decent. Despite the energies predicted by the model of García-Sánchez et al. also appears to be proportional to DFT energies, it shows a notably more repulsive behavior than that by DFT. In contrast, the energies predicted by the García-Pérez et al. have a relatively weaker correlation with the DFT energies. For instance, while the force field predicts an energy of approximately -2200 K, the DFT energies can range from -300 K to -3500 K. Although a conclusive evaluation with respect to the accuracy of these force fields still remains unknown, D2FF was identified to reasonably describe the atomic-level energy landscape of CO2 adsorbed in zeolite CHA. This outcome is in fact not surprising given D2FF was developed originally based upon DFT calculations.39 Recently, the development of force fields based on ab initio methods has drawn considerable attention, which has been demonstrated to lead to potential (including D2FF) that can accurately reproduce experimental properties.59–61 While DFT-derived force fields may seem promising, a proper choice of van der Waals functional is crucial and can be sometimes difficult. We therefore propose that a combination of DFT-computed reference energies with experimental measurements of adsorption properties including uptakes at a high pressure and diffusion coefficient may yield an accurate, physically sound, and transferrable force field.

2x103 García-Pérez et al. García-Sánchez et al. D2FF

1x103

DFT-calculated energy (K)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

0

-1x103

-2x103

-3x103

-4x103 -4x103

-3x103

-2x103

-1x103

0

1x103

2x103

Force field-computed energy (K)

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 30

Figure 8. Comparison of interaction energies of CO2 in zeolite CHA predicted by force fields of García-Pérez et al., García-Sánchez et al., and D2FF with that calculated by DFT with the nonlocal vdW-DF2 functional. A parity line is provided for the guide to the eyes.

4. Conclusion In this study, we have demonstrated that the predictions of adsorption properties made by available force fields, which have been widely adopted under a transferable assumption in screening studies, can result in a large discrepancy of more than two orders of magnitude. Such a discrepancy was also found to result in a non-negligible uncertainty in the evaluation of materials’ performance in CCS, evidently suggesting the importance of selecting suitable force fields in future large-scale studies. It was found that the difference in the repulsive characteristics of force fields (i.e., primarily the repulsive behavior of the van der Waals interactions) led to a large variation in predictions, in particular for geometrically confined zeolites with Di ≤ 5 Å. In contrast, different force fields generally yield a consistent prediction in Henry coefficients when structures are less confined. Interestingly, those structures that were included in the potential parameterization and/or validation were all fallen within the latter region. While including a diverse set of structures in the potential development may be less feasible, we have proposed several other properties that may help facilitate probing the repulsive nature of intermolecular interactions. Saturation loading (or uptakes at a high pressure) was identified to be an important property but has been largely overlooked. Such property indicates the pore volume available to gas molecules should be inversely correlated to the repulsive characteristics (i.e., more repulsive, less available volume). Mixture isotherms can also be used to validate the developed potential although they may be rather indirect. Moreover, gaseous diffusion coefficients may be another critical reference to be included, which are expected to be largely sensitive to confined space. Additionally, to ensure that physically sound force field is achieved for capturing atomic-level energy-landscape of adsorption, ab initio calculations such as DFT to compute interaction energies as references for force field development are suggested. Overall, the outcomes of this study are instrumental to the future development of force fields. Having accurate and transferable potentials can allow accurate identification of promising candidates from large-scale screening studies, thus pushing forward the discovery of new materials toward more energy-efficient and cost-effective applications. It should be noted that, although this study have only focused on all-silica zeolites, suggestions made from this study should be also generally applicable to other classes of nanoporous materials.

ACS Paragon Plus Environment

22

Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Author Contributions The study was developed and completed through contributions by all authors.

Supporting Information Available The Supporting Information includes additional tables and figures referred in the main article. Additional tables are shown for force field parameters and Spearman rank coefficients for uptakes predicted by different force fields. Additional figures are also shown for full-scale plots of Henry coefficient ratio against Di, plots for absolute Henry coefficients, Henry coefficient ratio plots w.r.t that predicted by the model of García-Pérez et al., and Henry coefficient ratio plots against structural features (i.e. density and accessible volume ratio). This material is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgement The authors thank the Ohio Supercomputer Center (OSC)62 for providing computational resources to carry out this work. J.R.Lim, C.-T.Yang, and L.-C.Lin also gratefully acknowledge the financial support from the Ohio State University.

ACS Paragon Plus Environment

23

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 30

References (1)

Solomon, S.; Plattner, G.-K.; Knutti, R.; Friedlingstein, P. Irreversible Climate Change due to Carbon Dioxide Emissions. Proc. Natl. Acad. Sci. 2009, 106, 1704–1709.

(2)

Hansen, A. J.; Johnson, D.; Lacis, A.; Lebedeff, S.; Lee, P.; Rind, D.; Russell, G. Climate Impact of Increasing Atmospheric Carbon Dioxide. Science 1981, 213, 957–966.

(3)

Chu, S. Carbon Capture and Sequestration. Science 2009, 325, 1599.

(4)

Ramezan, M.; Skone, T. J.; Nsakala, N.; Lilijedahl, G. Carbon Dioxide Capture from Existing Coal-Fired Power Plants; DOE/NETL-401/110907; 2007.

(5)

Coninck, H. C.; Loos, M. A.; Metz, B.; Davidson, O.; Meyer, L. IPCC Special Report on Carbon Dioxide Capture and Storage; Cambridge: New York, 2005.

(6)

Rochelle, G. T. Amine Scrubbing for CO2 Capture. Science 2009, 325, 1652–1654.

(7)

Wilmer, C. E.; Farha, O. K.; Bae, Y.-S.; Hupp, J. T.; Snurr, R. Q. Structure–property Relationships of Porous Materials for Carbon Dioxide Separation and Capture. Energy Environ. Sci. 2012, 5, 9849.

(8)

Yazaydin, A. O.; Snurr, R. Q.; Park, T.-H.; Koh, K.; Liu, J.; Levan, M. D.; Benin, A. I.; Jakubczak, P.; Lanuza, M.; Galloway, D. B.; Low, J. J.; Willis, R. R. Screening of Metal Organic Frameworks for Carbon Dioxide Capture from Flue Gas Using a Combined Experimental and Modeling Approach. J. Am. Chem. Soc. 2009, 131,18198–18199.

(9)

D’Alessandro, D. M.; Smit, B.; Long, J. R. Carbon Dioxide Capture: Prospects for New Materials. Angew. Chemie - Int. Ed. 2010, 49, 6058–6082.

(10)

Bhown, A. S.; Freeman, B. C. Analysis and Status of Post-Combustion Carbon Dioxide Capture Technologies. Environ. Sci. Technol. 2011, 45, 8624–8632.

(11)

Huck, J. M.; Lin, L.-C.; Berger, A. H.; Shahrak, M. N.; Martin, R. L.; Bhown, A. S.; Haranczyk, M.; Reuter, K.; Smit, B. Evaluating Different Classes of Porous Materials for Carbon Capture. Energy Environ. Sci. 2014, 7, 4132–4146.

(12)

Lin, L.-C.; Berger, A. H.; Martin, R. L.; Kim, J.; Swisher, J. A.; Jariwala, K.; Rycroft, C. H.; Bhown, A. S.; Deem, M. W.; Haranczyk, M.; et al. In Silico Screening of Carbon-Capture Materials. Nat. Mater. 2012, 11, 633.

ACS Paragon Plus Environment

24

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(13)

Flanigen, E. M.; Bennett, J. M.; Grose, R. W.; Cohen, J. P.; Patton, R. L.; Kirchner, R. M.; Smith, J. V. Silicalite, a New Hydrophobic Crystalline Silica Molecular Sieve. Nature

1978, 271, 512–516. (14)

Cruciani, G. Zeolites upon Heating: Factors Governing Their Thermal Stability and Structural Changes. J. Phys. Chem. Solids 2006, 67, 1973–1994.

(15)

Deem, M. W.; Pophale, R.; Cheeseman, P. A.; Earl, D. J. Computational Discovery of New Zeolite-like Materials. J. Phys. Chem. C 2009, 113, 21353–21360.

(16)

Pophale, R.; Cheeseman, P. A.; Deem, M. W. A Database of New Zeolite-like Materials. Phys. Chem. Chem. Phys. 2011, 13, 12407.

(17)

International Zeolite Association (IZA) databaes. http://www.iza-structure.org/databases/. (accessed March 5, 2018)

(18) Jeong, W.; Kim, J. Understanding the Mechanisms of CO2 Adsorption Enhancement in Pure Silica Zeolites under Humid Conditions. J. Phys. Chem. C 2016, 120, 23500–23510. (19)

Hasan, M. M. F.; First, E. L.; Floudas, C. A. Cost-Effective CO2 Capture Based on in Silico Screening of Zeolites and Process Optimization. Phys. Chem. Chem. Phys. 2013, 15, 17601.

(20)

Martin, R. L.; Willems, T. F.; Lin, L.-C.; Kim, J.; Swisher, J. A.; Smit, B.; Haranczyk, M. Similarity-Driven Discovery of Zeolite Materials for Adsorption-Based Separations. ChemPhysChem 2012, 13, 3595–3597.

(21)

First, E. L.; Gounaris, C. E.; Floudas, C. A. Predictive Framework for Shape-Selective Separations in Three-Dimensional Zeolites and Metal-Organic Frameworks. Langmuir

2013, 29, 5599–5608. (22)

Curtarolo, S.; Hart, G. L. W.; Nardelli, M. B.; Mingo, N.; Sanvito, S.; Levy, O. The HighThroughput Highway to Computational Materials Design. Nat. Mater. 2013, 12, 191–201.

(23)

Kim, J.; Lin, L.-C.; Martin, R. L.; Swisher, J. A.; Haranczyk, M.; Smit, B. Large-Scale Computational Screening of Zeolites for Ethane/ethene Separation. Langmuir 2012, 28, 11914–11919.

(24)

Nugent, P.; Giannopoulou, E. G.; Burd, S. D.; Elemento, O.; Giannopoulou, E. G.; Forrest, K.; Pham, T.; Ma, S.; Space, B.; Wojtas, L.; Eddaoudi, M.; Zaworotko, M. J..

ACS Paragon Plus Environment

25

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 30

Porous Materials with Optimal Adsorption Thermodynamics and Kinetics for CO2 separation. Nature 2013, 495, 80–84. (25)

Tagliabue, M.; Farrusseng, D.; Valencia, S.; Aguado, S.; Ravon, U.; Rizzo, C.; Corma, A.; Mirodatos, C. Natural Gas Treating by Selective Adsorption: Material Science and Chemical Engineering Interplay. Chem. Eng. J. 2009, 155, 553–566.

(26)

Peng, X.; Cao, D. Computational Screening of Porous Carbons, Zeolites, andMetal Organic Frameworks for Desulfurization andDecarburizat Ion of Biogas, Natural Gas, and Flue Gas. AIChE J. 2013, 59, 2928-2942.

(27)

First, E. L.; Hasan, M. M. F.; Floudas, C. A. Discovery of Novel Zeolites for Natural Gas Purification Through Combined Material Screening and Process Optimization. AIChE J.

2014, 60, 1767-1785. (28)

Braun, E.; Zurhelle, A. F.; Thijssen, W.; Schnell, S. K.; Lin, L.-C.; Kim, J.; Thompson, J. A.; Smit, B. High-Throughput Computational Screening of Nanoporous Adsorbents for CO 2 Capture from Natural Gas. Mol. Syst. Des. Eng. 2016, 1, 175–188.

(29)

Kim, J.; Abouelnasr, M.; Lin, L.-C.; Smit, B. Large-Scale Screening of Zeolite Structures for CO2 Membrane Separations. J. Am. Chem. Soc. 2013, 135, 7545–7552.

(30)

Kim, J.; Maiti, A.; Lin, L.-C.; Stolaroff, J. K.; Smit, B.; Aines, R. D. New Materials for Methane Capture from Dilute and Medium-Concentration Sources. Nat. Commun. 2013, 4, 1694–1697.

(31)

Sikora, B. J.; Wilmer, C. E.; Greenfield, M. L.; Snurr, R. Q. Thermodynamic Analysis of Xe/Kr Selectivity in over 137 000 Hypothetical Metal–organic Frameworks. Chem. Sci.

2012, 3, 2217-2223. (32)

Ryan, P.; Farha, O. K.; Broadbelt, L. J.; Snurr, R. Q. Computational Screening of MetalOrganic Frameworks for Xenon/Krypton Separation. AIChE J. 2011, 57, 1759-1766.

(33)

Jameson, C. J.; Jameson, A. K.; Lim, H.-M. Competitive Adsorption of Xenon and Krypton in Zeolite NaA: 129Xe Nuclear Magnetic Resonance Studies and Grand Canonical Monte Carlo Simulations. J. Chem. Phys. 1997, 107, 4364–4372.

(34)

Vujic, B.; Lyubartsev, A. P. Transferable Force-Field for Modelling of CO2, N2, O2 and Ar in All Silica and Na+ Exchanged Zeolites. Model. Simul. Mater. Sci. Eng. 2016, 24, 26.

ACS Paragon Plus Environment

26

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(35)

Bai, P.; Tsapatsis, M.; Siepmann, J. I. TraPPE-Zeo: Transferable Potentials for Phase Equilibria Force Field for All-Silica Zeolites. J. Phys. Chem. C 2013, 117, 24375-24387.

(36)

García-Pérez, E.; Parra, J. B.; Ania, C. O.; García-Sánchez, A.; Van Baten, J. M.; Krishna, R.; Dubbeldam, D.; Calero, S. A Computational Study of CO2, N2, and CH4 Adsorption in Zeolites. Adsorption 2007, 13, 469-476..

(37)

Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field. J. Phys. Chem. B 2004, 108, 1255-1266.

(38)

García-Sánchez, A.; Ania, C. O.; Parra, J. B.; Dubbeldam, D.; Vlugt, T. J. H.; Krishna, R.; Calero, S. Transferable Force Field for Carbon Dioxide Adsorption in Zeolites. J. Phys. Chem. C 2009, 113, 8814-8820.

(39)

Fang, H.; Kamakoti, P.; Zang, J.; Cundy, S.; Paur, C.; Ravikovitch, P. I.; Sholl, D. S. Prediction of CO 2 Adsorption Properties in Zeolites Using Force Fields Derived from Periodic Dispersion-Corrected DFT Calculations. J. Phys. Chem. C 2012, 116, 1069210701.

(40)

Dubbeldam, D.; Torres-Knoop, A.; Walton, K. S. On the Inner Workings of Monte Carlo Codes. Mol. Simul. 2013, 39, 1253-1292.

(41)

Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q. RASPA: Molecular Simulation Software for Adsorption and Diffusion in Flexible Nanoporous Materials. Mol. Simul.

2016, 42, 81–101. (42) Harris, J. G.; Yung, K. H. Carbon Dioxide’s Liquid-Vapor Coexistence Curve and Critical Properties as Predicted by a Simple Molecular Model. J. Phys. Chem. 1995, 99, 12021– 12024. (43)

Smit, B.; Maesen, T. L. M. Molecular Simulations of Zeolites: Adsorption, Diffusion, and Shape Selectivity. Chem. Rev. 2008, 108, 4125–4184.

(44)

Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Elsevier: London, 2001.

(45)

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593.

ACS Paragon Plus Environment

27

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(46)

Page 28 of 30

Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Nardelli, M. B.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; et al. Advanced Capabilities for Materials Modelling with Q Uantum ESPRESSO. J. Phys. Condens. Matter 2017, 29, 465901.

(47)

Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. A Higher-Accuracy van Der Waals Density Functional. Phys. Rev. B 2010, 82, 081101.

(48)

Troullier, N.; Martins, J. A Straightforward Method for Generating Soft Transferable Pseudopotentials. Solid State Commun. 1990, 74, 613–616.

(49)

Potoff, J. J.; Siepmann, J. I. Vapor-Liquid Equilibria of Mixtures Containing Alkanes, Carbon Dioxide, and Nitrogen. AIChE J. 2001, 47, 7.

(50)

Zimmermann, N. E. R.; Haranczyk, M.; Sharma, M.; Liu, B.; Smit, B.; Keil, F. J. Adsorption and Diffusion in Zeolites: The Pitfall of Isotypic Crystal Structures. Mol. Simul.

2011, 37, 986–989. (51)

Spearman, C. The Proof and Measurement of Association between Two Things. Am. J. Psychol. 1904, 15, 72-101.

(52)

Myers, J. L.; Well, A.; Lorch, R. F. Research Design and Statistical Analysis; Routledge: New York, 2010.

(53)

Bhatia, S. K.; Myers, A. L. Optimum Conditions for Adsorptive Storage. Langmuir 2006, 22, 1688–1700.

(54)

Simon, C. M.; Kim, J.; Lin, L.-C.; Martin, R. L.; Haranczyk, M.; Smit, B. Optimizing Nanoporous Materials for Gas Storage. Phys. Chem. Chem. Phys. 2014, 16, 5499.

(55)

Martin, R. L.; Smit, B.; Haranczyk, M. Addressing Challenges of Identifying Geometrically Diverse Sets of Crystalline Porous Materials. J. Chem. Inf. Model. 2012, 52, 308–318.

(56)

Myers, A. L.; Prausnitz, J. M. Thermodynamics of Mixed‐gas Adsorption. AIChE J. 1965, 11, 121–127.

(57)

Rao, M. B.; Sircar, S. Thermodynamic Consistency for Binary Gas Adsorption Equilibria. Langmuir 1999, 15, 7258–7267.

(58)

Swisher, J. A.; Lin, L.-C.; Kim, J.; Smit, B. Evaluating Mixture Adsorption Models Using Molecular Simulation. AIChE J. 2013, 59, 3054–3064.

ACS Paragon Plus Environment

28

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(59)

Lin, L.-C.; Lee, K.; Gagliardi, L.; Neaton, J. B.; Smit, B. Force-Field Development from Electronic Structure Calculations with Periodic Boundary Conditions: Applications to Gaseous Adsorption and Transport in Metal-Organic Frameworks. J. Chem. Theory Comput. 2014, 10, 1477–1488.

(60)

Dzubak, A. L.; Lin, L.-C.; Kim, J.; Swisher, J. A.; Poloni, R.; Maximoff, S. N.; Smit, B.; Gagliardi, L. Ab Initio Carbon Capture in Open-Site Metal-Organic Frameworks. Nat. Chem. 2012, 4, 810–816.

(61)

Fang, H.; Demir, H.; Kamakoti, P.; Sholl, D. S. Recent Developments in First-Principles Force Fields for Molecules in Nanoporous Materials. J. Mater. Chem. A 2014, 2, 274– 291.

(62)

Ohio Supercomputer-Center. Columbus OH. 1987. http://osc.edu/ark:/19495/f5s1ph73. (accessed March 5, 2018)

ACS Paragon Plus Environment

29

The Journal of Physical Chemistry

TOC Graphic 10.0

Differences attributed primarily to the prediction in repulsion among force fields

Henry coefficient ratio (w.r.t. D2FF)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 30

5.0

MFI CHA

1.5 1.0

0.5 Vujic et al. Bai et al. García-Pérez et al. CLAYFF García-Sánchez et al. 0.1

IZA structures

ACS Paragon Plus Environment

30