A process-level perspective of the impact of molecular force fields on the computational screening of MOFs for carbon capture

Abstract

The question we pose in this study is to what extent the ranking of metal organic frameworks (MOFs) for adsorption-based carbon capture, and the selection of top performers identified in Pressure Swing Adsorption (PSA) process modelling, depends on the choice of the commonly available forcefields. To answer this question, we first generated distributions of CO₂ and N₂ adsorption isotherms via molecular simulation in 690 MOFs using six typical forcefields: the UFF or Dreiding sets of Lennard-Jones parameters, in combination with partial charges derived from ab initio calculations or by charge equilibration schemes. We then conducted a systematic uncertainty quantification study using PSA process-level modelling. We observe that: (i) the ranking of MOFs significantly depends on the choice of forcefield; (ii) partial charge assignment is the prevailing source of uncertainty, and that charge equilibration schemes produce results which are inconsistent with ab initio-derived charges; (iii) the choice of Lennard-Jones parameters is a considerable source of uncertainty. Our work highlights that is not really possible to obtain material rankings with high resolution using a single molecular modelling approach and that, as a minimum, some uncertainty should be estimated for the performance of MOFs shortlisted using high throughput computational screening workflows. Future prospects for computational screening studies are also discussed.

---

Broader context

To avoid the most severe impacts of climate change, we need to deploy carbon capture and sequestration (CCS) technologies on a scale that matches human industrial activity's emissions. However, the differences in scale, composition, and technical requirements between various emission sectors require technologies that can be optimised for various applications. Adsorption-based processes using metal organic frameworks (MOFs) for carbon capture offer a significant advantage because they can be tailored for specific CO₂ capture sources by optimising the adsorption cycle and MOF properties. While high throughput computational screening (HTCS) studies are routinely conducted to identify the best MOFs for various CCS applications, very few of these materials are tested at the lab-scale. In this article, we argue that a lack of accurate, reproducible, and consistent implementation of HTCS workflows is one of the primary technological barriers between theoretical and experimental studies. We demonstrate that the molecular parameter sets used to simulate CO₂ and N₂ adsorption significantly impact the final recommendations made using process modelling. Given the uncertainties we observe, the community must take steps to improve the current simulation strategy. We offer several recommendations in the article on how to enhance the reliability of predictions derived from HTCS studies.

Introduction

High throughput computational screening (HTCS) has emerged as an important strategy in the efforts to identify the most promising metal organic frameworks (MOFs) for carbon capture applications.^1,2 Consider two academic groups as an example: both are seeking to determine the best MOF for carbon capture from flue gas using computational screening. Their simulation protocols are identical, including the database of MOF candidates. To rank the performance of materials, both groups use pressure swing adsorption (PSA) process modelling, with an identical process configuration and parameters. They both operate under the assumption that the results obtained using their workflows can be used to inform experimentalists on which MOFs to synthesise and test. The only difference between the groups is the choice of force fields used to describe intermolecular interactions, which are required by molecular simulations to generate equilibrium adsorption data. In this article, we question whether the rankings produced by the two groups, and therefore the final recommendations of their screening studies, will be similar and to what extent.

Let us first look in more detail at the steps and the assumptions typically involved in HTCS workflows. Generally speaking, HTCS can be conducted at the material level or the process level. Material-level screening studies filter MOFs for a desired property. For example, the CO₂/N₂ adsorption selectivity, percent regenerability, and CO₂ working capacity are key performance indicators (KPIs) that have been used to organise MOF databases on several occasions.^2–4 However, for complex, dynamic processes such as PSA, it has been convincingly demonstrated that these properties alone do not provide enough information to determine the separation potential of an adsorbent.⁵ For that reason, process-level screening workflows must be utilised. This approach relies on atomistic Grand Canonical Monte Carlo (GCMC) simulations to obtain microscale properties of the adsorbent (in particular, equilibrium adsorption data). This data is then passed onto detailed PSA process models whereby MOFs are evaluated under realistic separation conditions.⁶ At this level, we are interested in KPIs which reflect the economic drivers of the separation, such as high CO₂ product purity and recovery, or low energy penalty (energy required to capture a unit amount of CO₂) and high productivity (how much CO₂ is captured per kg of adsorbent per second). These properties depend on the properties of the adsorbent as well as the variables of the PSA cycle (pressures of each step, duration of each step, etc.). However, there is no unique combination of the PSA cycle variables which satisfies all of the process-level KPIs simultaneously, and so the trade-off between competing objectives can be represented as a series of Pareto-optimal solutions known as a Pareto front. Comparing Pareto fronts for two different materials has therefore become an established protocol for ranking MOFs based on their performance.^7–11 Among the several design choices one is posed with when constructing a process-level simulation pipeline, the first (and perhaps most important) choice is how to generate the adsorption data using molecular simulations. To understand the importance of this step, it is instructive to briefly review the primary components of a typical GCMC simulation of adsorption.

The interactions between all atoms, and therefore the potential energy, must be accurately described in order to predict CO₂ and N₂ adsorption in MOFs. In classical GCMC, this is typically achieved by defining an appropriate set of equations and molecular parameters known as a forcefield (FF), whereby short-range repulsion/dispersion interactions are described by a Lennard-Jones (LJ) potential, and long-range electrostatic interactions by a coulombic potential.^12,13 Specialised FFs have been developed to accurately describe the LJ contributions in specific subclasses of MOFs,^12,14,15 however these lack the large-scale predictive capabilities needed for screening purposes. On the other hand, generic FFs such as Dreiding¹⁶ and the Universal Force Field (UFF)¹⁷ are transferable (albeit less accurate) as they describe the same types of framework atoms in many different materials with the same parameters. Therefore, Dreiding or UFF are the common and practical choices for HTCS of MOFs.^18–22 They are typically combined with the Transferable Potentials for Phase Equilibria (TraPPE)²³ FF (which describes the CO₂ and N₂ adsorbates) using mixing rules such as Lorentz–Berthelot. For the electrostatic interactions, partial charges must be assigned to each atom in the system. Several assignment schemes have been developed so far and have been reviewed elsewhere,²⁴ each differing in their level of accuracy, philosophy, and computational cost. The most accurate (but computationally expensive) approach would be to determine the electrostatic potential within a material using ab initio calculations and then assign point charges that best represents this electronic environment.²⁴ While many ab initio approaches exist,^25–33 density derived electrostatic and chemical (DDEC)^28–32 charges and Repeating Electrostatic Potential Extracted ATomic (REPEAT) charges³³ are considered to be the best options for periodic structures. Less accurate, but computationally inexpensive, charges can also be obtained using charge equilibration (Qeq)³⁴ schemes such as the Extended Qeq (EQeq)³⁵ and Periodic Qeq (PQeq)³⁶ methods, among others.^37,38 These approaches are based on the principle of electronegativity equalisation between bonded atoms, and have been used extensively in MOF screening studies despite their semiempirical nature.^18,39–42

Despite the consensus that ab initio charges are the preferred choice for HTCS, several practical considerations (such as the availability of computational resources) can motivate the choice of charge equilibration schemes. Additionally, it is currently nontrivial to obtain a set of LJ parameters capable of screening thousands of MOFs with high accuracy. Therefore, it is quite common to see an almost arbitrary combination of the LJ parameters and charge assignment schemes used in HTCS of MOFs. This raises some concerns on the quality of the FFs used for predicting equilibrium adsorption data in molecular simulations.^43–45 Several studies have endeavoured to understand how the choice of charge scheme^27,39,46,47 or LJ FF^44,48 changes the rankings of MOFs for CO₂ capture. However, these studies do not extend beyond a simple material-level analysis (i.e., up to step 2 in Fig. 1). In order to align with the advancements achieved in HTCS, the uncertainty embedded within the process-level KPIs must be understood (i.e., up to step 3 in Fig. 1), given that the adsorption data plays a pivotal role in PSA-based CO₂ capture processes.^5,49–51 Without this level of analysis, a shadow of doubt on the outcomes of theoretical studies will remain, which only serves to hinder the transition of the recommended materials to the experimental campaign and ultimately to practical, industrial systems.

Fig. 1 An approach to understanding the impact of molecular forcefields on the computational ranking of MOFs for CO₂ capture. In the first step, a subset of molecular forcefields (comprised of Lennard-Jones parameters and partial charge assignment schemes) are chosen for evaluation. In the second step, adsorption data is generated using different forcefields via molecular simulation. The rankings of MOFs are then compared using material-level metrics such as CO₂ uptake. Until now, the differences between forcefields have only been evaluated at this level. In the third step, as in the current study, the impact of one's choice in forcefield is evaluated using rigorous process simulation and optimisation.

To explore these issues, we carry out a multiscale HTCS where we assess the suitability of six FFs commonly chosen to simulate CO₂ and N₂ adsorption in MOFs. Our goal is to understand how the material rankings are impacted by this choice, not necessarily on how well the FFs reproduce experimental data, and so we deliberately avoid making any comparisons to experimental measurements of adsorption. As a process-level case study, we consider the separation of CO₂ (15 vol%) and N₂ (85 vol%) binary mixtures using a modified Skarstrom PSA cycle at conditions representative of dry flue gas streams from coal-fired power plants. We generate distributions of CO₂ and N₂ adsorption isotherms in 690 MOFs at multiple temperatures using LJ parameters from either Dreiding or UFF, in combination with partial charges assigned by the DDEC, EQeq, or Qeq schemes. We then determine the distributions in process-level performances for every MOF using a combination of high-fidelity PSA process modelling and machine learning. Using these data, we generate material ranking lists for each FF and reflect on the consistency of the results. Our analysis will show moderate-to-poor correlations between the process-level rankings obtained using different FFs, and that this disagreement stems primarily from the significant deviations which manifest at the material level of description. We then quantify the average uncertainties between FFs from their process-level responses and reveal that the prevailing source of uncertainty is the choice of charge scheme. We finally explore a potential pathway towards uncertainty mitigation in multiscale simulations of PSA-based CO₂ capture processes and discuss future prospects for HTCS of materials.

Methodology

Grand canonical Monte Carlo simulations

The computation-ready, experimental (CoRE) MOF 2014 DDEC database^52,53 is chosen as the starting point for adsorbent selection. To ensure a fair comparison in GCMC modelling approaches, only MOF structures that contain all the atom types present in both the UFF and Dreiding force fields are included. This filtering process results in a reduced dataset of 726 structures from the original 2932. A further 36 MOFs are removed due to the presence of unbound water molecules, counter-ions and/or hydrogen atoms with incorrect bond lengths/angles, leaving 690 computation-ready structures available for this study.

Simulations of CO₂ and N₂ adsorption are performed using the atomistic GCMC method implemented in RASPA.⁵⁴ For framework atoms, LJ parameters were taken from UFF¹⁷ or Dreiding.¹⁶ For CO₂ and N₂ adsorbates, the TraPPE FF is used.²³ For the coulombic interactions, three partial charge assignment schemes are considered here: DDEC,²⁸ EQeq,³⁵ and Qeq.³⁴ DDEC charges are taken from the original CoRE MOF 2014 DDEC database, while the EQeq and Qeq partial charges are calculated using the EQeq v1.1.0⁵⁵ and RASPA v2.0.47 simulation codes, respectively.

For each adsorbent, single-component isotherms are generated at 273 K, 298 K, and 323 K on a logarithmic pressure scale between 10⁻³–10 bar for CO₂ and 10⁻³–1 bar for N₂. As there are 6 possible FF combinations, this results in 18 isotherms for both CO₂ and N₂ using 378 unique adsorption points. In total, 24 840 isotherms are generated across 690 MOFs which span different pressures, temperatures, charge assignment methods, and generic forcefield parameters. We refer to the data generated using this subset of MOFs as the charge-dependent, reproducible, accessible, forcefield-dependent, and temperature-dependent exploratory database of isotherms, or CRAFTED isotherms for brevity. To provide an initial seed for future benchmarking efforts, we make the CRAFTED database (v1.1.1)⁵⁶ and scripts available to the broader scientific community in a dedicated Zenodo repository at https://zenodo.org/record/7689919. This data was obtained using a HTCS workflow that has been developed⁵⁶ and validated elsewhere.^57,58 For more technical details on the GCMC simulations, see the ESI,† supplementary note 1.

Equilibrium adsorption model

Process-level simulations require an analytical model to determine the equilibrium molar loadings of CO₂/N₂ mixtures over various pressure and temperature domains. The dual-site Langmuir (DSL) model (eqn (1)) is well suited for this purpose. Within the DSL approach, competitive adsorption may be calculated by a simple extension of the single-component DSL models;⁵⁹ many practical systems of interest are known to be described by this behaviour.^7,60,61Here, i refers to adsorbate (either CO₂ or N₂), j refers to adsorption sites 1 or 2, qⁱ_s,j [mmol g⁻¹] is the saturation capacity, p_i [bar] is the pressure of component i, qⁱ_0,j [bar⁻¹] is the pre-exponential factor, and ΔH^ads_i,j [J mol⁻¹] is the heat of adsorption. The DSL model parameters are determined by nonlinear least-squares regression against the GCMC adsorption data at all three temperatures using the lmfit.minimize() function in python.⁶² We apply the following constraints during the fitting to obtain physically meaningful DSL model parameters. We enforce that the saturation capacity of each site is equal for both CO₂ and N₂ to satisfy the thermodynamic consistency requirements.⁶³ We designate the CO₂ adsorption site 1 as being the stronger adsorbing site by imposing the condition. This constraint reflects the energetic heterogeneity of CO₂ adsorption. For N₂, in cases where the adsorption isotherm is linear over the pressure range, we adopt the Equal Energy Site (EES) formulation of the extended DSL model.^64,65 This formulation recognises the energetic homogeneity of N₂ adsorption by setting the interaction energy parameters (i.e., and ) to the same value at both sites.^51,66 In cases where N₂ adsorption deviates substantially from linearity, this constraint is relaxed, and the fitting protocol used to model CO₂ adsorption is utilised. This approach has been validated using binary mixture GCMC simulations (see supplementary note 1.3, ESI†).

PSA simulation and optimisation strategy

We consider a modified Skarstrom PSA cycle which consists of five steps: (1) pressurisation; (2) adsorption; (3) heavy reflux; (4) counter-current depressurisation; and (5) light reflux, as shown in Fig. 2. The cycle variables consist of high (adsorption) and low (evacuation) pressures (p_H and p_L [bar]), durations of the pressurisation, adsorption, and counter-current depressurisation steps (t_pres, t_ads, and t_cnCDepres [s]), velocity of the feed (v_feed [m s⁻¹]), and light and heavy reflux ratios (α_LR and α_HR [−]). Other important properties highlighted in Fig. 2 are the porosity of the bed and pellet (ε_B and ε_p [−]), and the pellet radius (r_p [m]). The role and purpose of individual steps are described in detail elsewhere.⁷ The main feature of this cycle is that it utilises a dual reflux configuration, which avoids the need to pull deep vacuum below the practical limits of industrial vacuum pumps.⁶⁷ This presents an important advantage over the more commonly used 4-step VSA cycle with light product pressurisation, which often requires unphysically low evacuation pressures to achieve high CO₂ purities and recoveries.⁵¹

Fig. 2 Schematic of the five-step modified Skarstrom cycle. From left to right, the pressurisation, adsorption, heavy reflux, counter-current depressurisation, and light reflux steps are shown. Some of the key cycle variables and adsorbent pellet properties are highlighted in red and are described in the main text.

To accurately simulate the separation of CO₂/N₂, partial differential equations which describe the mass, energy, and momentum transfer phenomena taking place in the adsorption column are required. A complete description of the governing equations, model parameters, boundary conditions and mathematical criterion used to determine cyclic steady state are provided in supplementary note 2 (ESI†), while other technical details of the PSA simulation are provided elsewhere.^7,68 We summarise here the key assumptions adopted in our process model as the following: (1) the Ideal Gas law describes the gas phase; (2) the gas phase and solid phase are in thermal equilibrium; (3) mass transfer is controlled by molecular diffusion in the macropores, which is valid for most microporous adsorbents with pore openings greater than 4 Å.⁶⁹ Furthermore, the mass transfer rate is described by the linear driving force (LDF) approximation; (4) an axially dispersed plug flow model is employed, meaning there are no radial distributions in concentration, pressure, and temperature for both the solid and gas phase; (5) the column is operated adiabatically; (6) the Ergun equation describes the axial pressure drop. These assumptions are largely consistent with the assumptions of several previous studies.^7–11

Four key performance indicators (KPIs) are extracted from the PSA simulations, namely: CO₂ purity [−], CO₂ recovery [−], energy penalty [kWh ton_CO2⁻¹], and productivity [mol_CO2 kg_ads⁻¹ s⁻¹]. Regulatory bodies such as the US Department of Energy have specified 95% purity and 90% recovery targets for CO₂ capture processes. Viable adsorbents must therefore be able to achieve this product specification in any given PSA cycle. For the purposes of our investigation, we relax the purity requirement to 90%,^7,66 and simply refer to the 90% purity and recovery targets as carbon capture and sequestration (CCS) constraints. The energy penalty corresponds to the amount of energy required to capture a unit of CO₂, whereas productivity specifies how much CO₂ a unit of the adsorption bed (e.g., a kg or m³) can capture per unit of time. These KPIs are related to the operating and capital costs of the separation, and have been widely adopted in previous screening studies as the indicative metrics of material performance.^7–11,70 Mathematical definitions for each KPI are provided in supplementary note 2 (ESI†).

The performance of the PSA cycle is optimised by coupling the PSA process model with the Non-dominated Sorting Genetic Algorithm (NSGA-II).^71–73 11 decision variables (DVs) are initially considered for the process optimisation. These included all of the PSA cycle variables as well as ε_B, ε_p and r_p. Parameters ε_B, ε_p and r_p are included as DVs in light of the recent work by Farmahini et al.⁷⁴ who demonstrated that significant performance gains can be obtained by optimising properties of the pellet as well as the cycle variables. We then reduced the dimensionality of the optimisation landscape to 7 DVs by fixing the values of p_L, t_pres, t_cnCDepres, and α_HR after some preliminary sensitivity analysis. This decision is motivated by the fact that these 4 parameters contributed very little to the overall performance of the cycle; they either converged to the boundaries defined by physical limitations, or they varied within a narrow window of the optimisable range (see supplementary note 3, ESI†). A summary of the final DVs considered in this work and their operating ranges are provided in Table 1.

Table 1 Lower and upper bounds for decision variables used in the optimisation of the modified Skarstrom cycle

	p H [bar]	t ads [s]	α LR [−]	v feed [m s^−1]	ε B [−]	ε p [−]	r p [m]
a Competitive adsorption of CO2/N2 mixtures are well described by the extended DSL model up to 5 bar, see Fig. S1 (ESI).
Lower bound	1	5	0.01	0.1	0.35	0.3	0.5 × 10^−3
Upper bound	5	500	0.2	2.0	0.45	0.7	2.5 × 10^−3

---

We conduct two separate optimisation stages in this study. The first is an unconstrained maximisation of the CO₂ purity and recovery. For this stage, the genetic algorithm is terminated after 70 generations. The second stage is a multi-objective optimisation, whereby we maximise the productivity while simultaneously minimising the energy penalty subject to the two CCS constraints: CO₂ purity ≥0.9 and CO₂ recovery ≥0.9. For this stage, the genetic algorithm is terminated after 250 generations. We confirm that 70 and 250 generations for the CO₂ purity-recovery and energy-productivity optimisations, respectively, are sufficient to guarantee convergence of the Pareto fronts.

Machine learning surrogate model for the prediction of PSA performance indicators.

We developed a machine learning (ML) surrogate model of the modified Skarstrom cycle which can predict the process KPIs for any arbitrary combination of DVs and material properties. We use this surrogate model to expedite the calculation of CO₂ purity-recovery and energy penalty-productivity Pareto fronts. Similar to, and guided by, our previous efforts in this domain,⁶⁶ we use an artificial neural network (ANN) as our choice of ML model. There are overall 7 DVs and 14 material-specific parameters (crystal density ρ_s [kg m⁻³], crystal specific heat capacity c_p,s [J kg⁻¹ K⁻¹], and the 12 DSL parameters which describe CO₂ and N₂ adsorption), resulting in 21 inputs, while the outputs of interest are the 4 KPIs (i.e., CO₂ purity, CO₂ recovery, energy penalty, and productivity). Appropriate training limits are defined for each of the material-specific parameters to sufficiently cover the distributions of CRAFTED MOF material properties, and training limits for the 7 DVs are defined as the optimisable ranges given in Table 1.

Training data is generated by a combination of Latin Hypercube sampling (LHS) of the input phase space as well as a guided search of high quality, Pareto-optimal points through a bootstrap optimisation technique. The data obtained using these sampling techniques is then cleaned and pre-processed by applying log-like transformations to each of the KPIs.⁷⁵ Next, the data is partitioned into train/test/validation datasets using a ratio of 90/5/5. MATLAB's machine-learning toolbox is subsequently used to train a dense feedforward ANN model architecture of 3 hidden layers with 45 neurons per hidden layer by Levenberg–Marquardt back-propagation. The ANN model predictive performance is evaluated by calculating the adjusted coefficient of determination (R_adj²) for all KPIs using the test dataset. Model training and refinement is terminated once predictions of both the individual cycle configurations and fully resolved Pareto fronts no longer improved (Fig. 3). Additional details on the surrogate model training protocol and on how the CRAFTED MOF material properties are calculated are provided in the supplementary note 4 (ESI†).

Fig. 3 The ANN model prediction accuracy compared to the high-fidelity PSA process model. Left subplots: Parity plots showing the ANN model predictions of individual cycle configurations for each process KPI, along with the corresponding R_adj² values. Each ANN model is trained on 300 000+ datapoints. Right subplot: Examples of energy-productivity Pareto front predictions using the ANN model (solid lines) compared to the Pareto fronts obtained from process simulations (filled circles) for 7 different hypothetical adsorbent materials (indicated by different colours). Hypothetical adsorbent materials are characterised by random combinations of the material-specific parameters. Shaded region indicates an error of ±5%.

Uncertainty metrics for high performing adsorbents

High performing adsorbents are characterised by their ability to minimise the cost of CO₂ capture while meeting the CCS requirements. It is not a priori evident what combination of material properties leads to these requirements being satisfied, so we conduct preliminary purity-recovery optimisations in two stages. Low performing adsorbents are first identified by their ANN-generated Pareto fronts and discarded from our subsequent analysis. We define a low performing material as one which does not satisfy the CCS constraints for any FF. For the 218 materials that remain after this pre-screening, high quality Pareto data is generated by refining the ANN results with the PSA process model for an additional 10 NSGA-II generations. This mitigates any potential bias in the ML process and ensures that the conclusions drawn are established using accurate process simulation data. Moving forwards, we only consider the 218 MOFs in our uncertainty calculations. We designate this subset of MOFs as CRAFTED-u.

We seek to quantify the level of agreement between FFs at the molecular and process simulation levels. To do so, we make use of a number of statistical metrics which are both local, i.e., relating to a single material, and global, i.e., relating to the entire CRAFTED-u MOF database. Let us first focus on molecular simulations, used to generate equilibrium adsorption data. In addition to the comparison of individual adsorption isotherms generated from different FFs, we also employ the Spearman correlation coefficient, ρ, to measure the global extent of the correlation between adsorbate uptakes predicted by different FFs. Perfect positive correlations result in ρ = 1, perfect negative correlations in ρ = −1, and no correlation in ρ = 0.

At the process level, for each CRAFTED-u MOF, the uncertainty between Pareto fronts that emerges from the use of different molecular FFs to generate the adsorption data as input, is quantified using the hypervolume. The hypervolume, ξ, measures the area enclosed by all solutions on the Pareto front and a user-defined reference point r.^76,77 We measure ξ by querying N = 10⁴ uniformly distributed random points within the Pareto phase space bounded by r, where r is determined such that the entire Pareto front for every FF is covered in the sampling. The difference in hypervolumes between FFs i and j, which we call the hypervolume error Δξ^j_i, is taken as the uncertainty measure for a single material. The process of calculating Δξ^j_i is visually represented in Fig. 4.

Fig. 4 High throughput computational screening workflow with uncertainty quantification. For a single CRAFTED MOF, comparing two different molecular forcefields i and j involves the following steps. In step (1) the CO₂ and N₂ adsorption isotherms are predicted using molecular simulations. In step (2), the GCMC adsorption data is used to generate CO₂ purity-recovery Pareto fronts or energy penalty-productivity Pareto fronts using the ANN surrogate model. If the CCS constraints are satisfied for at least one forcefield, then all of the ANN-generated Pareto fronts are refined using the PSA process model. If the CCS constraints are not satisfied for any forcefield combination, then the MOF is discarded from our uncertainty quantification study. In step (3), the hypervolume ξ of each Pareto front i and j is determined by stochastically sampling the Pareto phase space with 10⁴ points, i.e., the blue and red shaded regions in the right-most subplots. The hypervolume error between forcefields i and j for a single material at the Pareto level, Δξⁱ_j, is taken as the difference in these two shaded regions, as shown by the grey shaded region in the right-bottom subplot.

Now, in order to quantify the global uncertainties between two different FFs i and j across all CRAFTED-u MOFs, we must determine a single numerical value from the full distribution of Δξ^j_i values. We take the mean of the hypervolume error distribution, , as the indicative metric of global uncertainty. Qualitatively we would also like to understand if FF i has a propensity to generate adsorption behaviours that produce better process performances overall compared to FF j. Establishing this connection between the adsorption behaviours and process-level responses helps to reveal the biases in performances that may arise in particular FFs. We therefore introduce a normalised metric, to understand the general process-level behaviours of different FFs across many materials which may dominate different areas of KPI phase space. The protocol for calculating for FF i is as follows. For a single MOF, six Pareto fronts are generated for each of the six possible FFs. One of these FFs generates a Pareto front which performs the best, and therefore has the greatest hypervolume. We treat this hypervolume as the upper limit of performance (i.e., ξ^max = 100%) and normalise the hypervolume of FF i relative to ξ^max such that Ξ_i = ξ_i/ξ^max × 100. This operation is applied to every CRAFTED-u MOF, resulting in a distribution of Ξ_i values. Similar to the hypervolume error, we measure the central tendency of the Ξ_i distribution by a single numerical value. In this case, we take the median, Ξ_i. If , then FF i tends to dominate over FF j overall. Note that the mean of Δξ^j_i and the median of Ξ_i were determined to be the most appropriate statistical measures of central tendency given the spread and skewness of the underlying data distributions. They are therefore the best metrics to indicate the overall trends of the CRAFTED-u MOF dataset.

Finally, we measure the correlation between process-level rankings using the Spearman coefficient. Note that, as a matter of convention, we indicate that any of the metrics discussed above are being calculated with respect to a particular FF by the subscript notation and are being compared to another FF using the superscript notation. As an example, ξ_i refers to the hypervolume calculated for FF i, and ρ^j_i refers to the correlation in rankings obtained using FFs i and j.

Results & discussion

We begin our discussion by delving into the most crucial and practical consideration of this work, i.e., to what extent the selection of a particular FF impacts the process-level rankings of materials. We then consider specific adsorption patterns and the level of agreement between different FFs in predicting the uptake of CO₂ and N₂. The objective here will be to reason some of the important differences in molecular modelling approaches, seeking to understand the nature of the (dis)agreement by exploring correlations in simulated CO₂ and N₂ uptakes. Next, we quantify the uncertainty between different FFs at the process-level. For this, we will utilise a number of figures and statistical measures to map process-level responses onto the underlying adsorption patterns. Finally, we explore potential pathways towards uncertainty mitigation in multiscale simulations of PSA-based CO₂ capture processes and discuss future prospects for HTCS of materials.

Level of agreement between process-level rankings obtained using different molecular forcefields

In this section, we inquire whether the ranking of porous materials and the identification of top performers through process modelling depends on the use of a specific FF. For this, we must first generate a material ranking list for each molecular FF using some process-level performance metric. The energy-productivity Pareto front, subject to the CCS constraints, is often used to rank materials in multiscale HTCS studies.^7,66,70,74 However, the CCS constraints are not satisfied by every CRAFTED-u MOF across all FFs, and so it is not possible to evaluate every material behaviour using energy-productivity KPIs. We therefore ranked materials according to the maximum CO₂ purity they can achieve in the modified Skarstrom cycle (Fig. 2), subject to the constraint of CO₂ recovery ≥0.9, using the HTCS workflow shown in Fig. 4.

Then, we extracted the top 50 MOFs from each ranking, consistent with the philosophy of previous screening studies which typically seek to narrow the candidate list to only the top performers.^{2,19,21,22,70} We find that the material rankings are poorly-to-moderately positively correlated, and that the identity of the top performing MOFs can change significantly, depending on the FF. To demonstrate this point, let us make comparisons amongst the different FFs by choosing UFF with DDEC charges as the reference. Note that by doing so, we do not suggest that UFF + DDEC is more accurate in reproducing experimentally observed adsorption behaviours, it is simply a matter of convention. With this in mind, we identify the number of MOFs from each FF-specific ranking which overlap with the top 50 materials from the UFF + DDEC baseline and visualise the results in Fig. 5.

Fig. 5 The number of MOFs in the list of top 50 performing materials shared between the UFF + DDEC forcefield and other molecular forcefields. Materials are ranked by the maximum CO₂ purity value they can achieve in the modified Skarstrom PSA cycle, subject to CO₂ recovery ≥0.9. The Spearman coefficient, ρ, measures the strength of the correlation between rankings obtained using the UFF + DDEC forcefield and other molecular forcefields, and is reported above each bar plot.

Let us first focus on the comparison between UFF and Dreiding in Fig. 5, with each of these LJ parameter sets combined with the DDEC charges. We see that 35 MOFs are shared in the list of top performers between these two variants of the forcefields. Turning now to the combinations of UFF with other charge schemes, the greatest level of agreement is between DDEC and Qeq, with 36 MOFs appearing in both lists of the top 50 materials. It is also quite apparent that the EQeq charges provide a much poorer level of agreement, with only 27 MOFs appearing in the ranking obtained using the DDEC charges. We find that even if we take a different molecular FF as the baseline (Fig. S9, ESI†) or focus on a shorter list of candidate materials, i.e., the top 20 (Fig. S10, ESI†), that generally: (1) there is consistently better agreement in the rankings between Qeq and DDEC charges than in EQeq and DDEC charges; (2) depending on one's choice of the charge scheme, the overlap in the materials identified by any two different schemes ranges between 50% and 70% (Fig. S9, ESI†), not to mention the order of materials which changes significantly depending on the set of parameters; (3) the screening results remain sensitive to the choice of LJ parameters and the two lists of top materials produced by two different FFs may share approximately 70% of the materials. These tendencies are supported by the Spearman ranking correlation coefficients, which are provided above each bar plot in Fig. 5.

While Fig. 5 is a useful tool to visualise how the identity of the top performers can change with molecular FF, we do not yet understand how close the materials are in terms of their numerical performance. The two possible explanations for poor-to-moderate agreement between material rankings is either: (1) many materials perform similarly, and small differences in the performance estimates give rise to large differences in the relative rankings across different FFs, or (2) the process performances, and therefore the rankings, are highly variable parameters. In Fig. 6 we provide a snapshot of how the numerical performance values used to rank the top 10 materials for UFF + DDEC change with molecular FF. From this plot we identify some materials, such as QOVWOO and PEKKAS, which show fairly minor deviations using different FFs, and so are consistently ranked in the top 10 candidates. Other materials, such as VUNCAJ, always remain within the top 50 threshold. Finally, we have some materials with such a strong volatility in their process KPIs that they may fall far below the top 50 threshold, depending on the choice of FF. These constitute the largest proportion of materials, and in many instances the difference in performance is sufficient to classify the same material as suitable (in other words, meeting the CCS constraints) by one FF and not suitable by another FF. Note that these trends persist when the other FFs are taken as the baseline for comparison (Fig. S11, ESI†).

Fig. 6 Visualisation of how the process-level performance values for the top 10 MOFs ranked by UFF + DDEC change using other molecular forcefields. Materials are ranked by the maximum CO₂ purity they can achieve in the modified Skarstrom cycle, subject to the recovery of CO₂ ≥ 0.9. CSD reference codes for the top 10 MOFs ranked by the UFF + DDEC forcefield (QOVWOO is ranked 1st, and KAFYAT is ranked 10th) are shown in the legend. The dark grey shaded region indicates the process-level performance value which differentiates the top 50 materials in each forcefield. The light grey shaded region indicates the CCS constraint of CO₂ purity ≥0.9.

So far it appears that the process performances and material rankings depend on both the definition of the LJ and the electrostatic potentials, however it seems that the choice of charge scheme has a stronger impact (Fig. 5). In practice, (E)Qeq may be chosen over DDEC in HTCS^78–80 due to restricted computational resources, or they may be used to pre-screen large databases and shortlist top performers for further refinement using higher accuracy charge assignment schemes.^81,82 In an ideal scenario, DDEC charges are available or can be calculated for all candidate materials. In this case, the uncertainty embedded in one's choice of partial charge assignment would be eliminated as ab initio charges are the preferred level of molecular theory for HTCS of materials for carbon capture.² If we confine the scope of our analysis to this ideal scenario, the only additional degree of freedom remaining is how to model the LJ interactions. In Fig. 7(a), a direct comparison of the maximum CO₂ purity predicted by either UFF or Dreiding in combination with DDEC charges reflects this ideal case study. Apart from FUFREE, the MOFs in this top performing subset show similar process-level responses. While it is tempting to attribute this behaviour to minor scattering in the adsorption data, we show later in our discussion of the process-level uncertainties that this may not strictly be the case as different adsorption patterns can lead to the same CO₂ purity-recovery performances. On the contrary, differences between LJ FFs become much more pronounced at the energy-productivity Pareto level. For example, in Fig. 7(b) we rank materials according to the maximum productivity they can achieve on their energy-productivity Pareto front (subject to the CCS constraints), and make a similar comparison between the top performers identified by UFF + DDEC and Dreiding + DDEC. It is evident from Fig. 7 that the predictions between LJ FFs are substantially more uncertain when the energy-productivity process KPIs are used. This observation is notable because the energy-productivity Pareto front more accurately reflects the separation potential of an adsorbent and, unlike the choice of charge scheme, the uncertainty that arises between different LJ parameter sets is unavoidable as the choice of UFF vs Dreiding is still somewhat arbitrary.

Fig. 7 Visualisation of how the process-level performance values for the top 10 MOFs ranked by UFF + DDEC change using Dreiding + DDEC. (a) Materials are ranked by the maximum CO₂ purity they can achieve in the modified Skarstrom cycle, subject to the recovery of CO₂ ≥ 0.9. (b) Materials are ranked by the maximum productivity they can achieve, subject to the constraints of CO₂ purity and recovery ≥0.9. PARHEW, ZNGLUD01, FUFREE, and PUPNAQ do not meet the CCS constraints using Dreiding + DDEC and so their productivity KPIs are set to zero. CSD reference codes for the top 10 MOFs in each subplot are shown in the legends.

What then are the practical implications of these observations for HTCS of materials for CO₂ capture? At this point it is important to recognise that the correlations between rankings are in some cases much greater than zero, i.e., between UFF + DDEC and Dreiding + DDEC. This suggests that HTCS retains its utility as a guided search platform even in the presence of considerable uncertainties. Indeed, Chung et al.⁸³ and Boyd et al.⁸² have demonstrated this to good effect as studies which have experimentally verified the structures predicted by HTCS workflows. However, in light of our results, the preliminary conclusion of this section is that it is very likely that two separate studies utilising different levels of molecular theory will still come to different conclusions on what the top performing material candidates are. The picture that is emerging from the process-level results cannot be complete without a complementary understanding of the adsorption behaviours. Our objective in the next section is to therefore understand the nature of the disagreement between FFs at the material level. This will provide us with the molecular insights needed to reason why such differences in process performances arise.

Analysis of the underlying adsorption patterns and the extent of consistency in CO₂ and N₂ uptakes predicted using different molecular forcefields

Fig. 8 provides a birds-eye-view of the extent of agreement in simulated uptakes of CO₂ and N₂ between different FFs. Fig. 8(a, b, e and f) show the Spearman correlations in simulated uptakes for all pressure points. Subplots (a) and (e) show correlations between LJ FFs using fixed charges for CO₂ and N₂, respectively, while subplots (b) and (f) consider fixed LJ FF parameters and different charge schemes. Parity plots of simulated uptakes for both CO₂ and N₂ relevant to each of these scenarios are provided in Fig. 8(c, d, g and h) at 0.001 bar, i.e., in the Henry's regime. We begin the discussion by comparing the behaviours of both adsorbates, after which we will explore the differences between molecular FFs. Note that the results presented here are qualitatively similar at higher pressures (see supplementary note 6, ESI†).

Fig. 8 Correlations and parity plots of simulated uptakes of CO₂ and N₂ for different molecular forcefields. Subplots (a) and (b) show the Spearman correlations in simulated uptakes of CO₂ for different LJ forcefields and charge schemes, respectively. Subplots (c) and (d) show the parity in simulated CO₂ uptakes at low pressure using different LJ forcefields and charge schemes, respectively. Subplots (e) through to (h) make the same comparisons but for N₂. Note that each point in subplots (c), (d), (g) and (h) represents a unique material and forcefield. Clusters I and II are highlighted in subplots (c) and (h) respectively to aid in the discussion.

From Fig. 8(a, b, e and f), both CO₂ and N₂ there is a greater degree of scattering, and therefore a lower correlation, between different FFs at lower pressures. We note that CO₂ simulated uptakes exhibit a higher degree of scattering than N₂ and seem more sensitive to the electrostatics than to the parameters of the LJ FF, as seen by the stronger correlations in Fig. 8(a) compared to Fig. 8(b); for N₂, the opposite trend is observed. The observations above are easy to rationalise. It is known that the most favourable adsorption sites on a MOF will be occupied preferentially at lower pressures.⁴⁴ The strength of the interactions between the adsorbing species and the framework is mediated by the partial charges assigned to framework atoms and the cross-interactions between the adsorbate and adsorbent (and hence the parameters of the LJ FF). Therefore, a higher degree of scattering in the uptakes is expected in the low pressure regime. CO₂ also exhibits more pronounced uncertainties than N₂ because: (1) it has more LJ sites available to interact with the MOF, and (2) it has a larger quadrupole moment which strongly interacts with the partial charges assigned to framework atoms.⁸⁴ Finally, as the electrostatic contributions in CO₂ adsorption are more significant than in N₂ adsorption, the variation in the charge schemes produces a more pronounced scattering in the results compared to variations in the LJ parameters.

Now, we explore in more detail individual adsorption behaviours that contribute to the overall picture in Fig. 8. Intuitively, we can expect three primary scenarios for when predictions between FFs deviate significantly from each other: (1) when adsorption occurs in MOFs with small pores; (2) when the LJ parameters for framework atoms are significantly different between UFF and Dreiding; and (3) when charge equilibration methods fail to reproduce DDEC charges. In addition, it is possible that there are combinations of these scenarios for particular MOFs leading to strongly coupled effects. The aim of the following discussion is to explain how these FF deviations arise, and in doing so to provide more insights into the molecular origins of the correlation results in Fig. 8(a, b, e and f). To achieve this, we explore each of these three scenarios, touching on both the general trends and some material specific deviations.

MOFs with smaller pores, and therefore a higher density of framework atoms, are more sensitive to the LJ parameters. These materials are typically associated with lower adsorbate uptakes, and so we generally see higher scattering in the data at lower loadings, as shown in e.g., Fig. 8(c). This is further demonstrated in Fig. S12 (ESI†), where the relative uptake of CO₂ at 0.1 bar is plotted against the accessible void fraction of MOFs. In specific cases, confinement effects may also lead to very strong deviations in the adsorptive behaviours depending on the choice of the LJ parameters. For instance, OSOMIT, HELCUY, and GIMVAA – materials which contain Zn and have very narrow pore spaces – represent such exceptional cases. As an example, let us consider the adsorption isotherms of GIMVAA, as shown in the top left panel of Fig. 9. The first feature to note about the behaviour of this material is the very strong uptake observed for the UFF + EQeq combination of parameters (green line). This is associated with the EQeq charges assigned to metal atoms in this structure, which are substantially higher than the charges assigned by the DDEC and Qeq schemes. Focusing on the other data in this panel (blue and red lines), we note that only UFF predicts an appreciable uptake of CO₂, and much lower uptake for any other combination of parameters. This behaviour stems primarily from the difference in LJ diameters assigned to Zn by UFF and Dreiding (σ^UFF_Zn = 2.462 Å and σ^Dre_Zn = 4.045 Å). Due to the larger size of the framework atoms, the cross-potential interactions between the adsorbate and adsorbent are unfavourable in the small pores of GIMVAA for Dreiding, and favourable for UFF. This leads to negligible CO₂ uptake when using Dreiding-based parameters and non-negligible uptake when the UFF models are employed. Returning to the outlier behaviour of the models based on EQeq charges, GIMVAA is representative of several materials in this class, forming a cluster of green points (cluster I) far above the parity line in Fig. 8(c).

Fig. 9 Distribution of CO₂ adsorption isotherms predicted by different molecular forcefields for GIMVAA, BEPNIV, AYUWUM, and EBOTOF. Error bars are shown at each simulated pressure point. MOF structures are visualised using iRASPA to the right of each isotherm subplot.

Let us now consider the second scenario, when the LJ parameters for framework atoms are significantly different between UFF and Dreiding. If a MOF contains atom types for which the depth of the well in the interaction potential, ε, between UFF and Dreiding is very different, we can anticipate some sensitivity of the adsorption behaviour to the choice of LJ FF.⁴⁸ From Fig. 8(c and g), we see that UFF generally predicts higher adsorbate uptakes. This can be tentatively attributed to the fact that many (approximately 65%) of the CRAFTED-u MOFs contain Zn as the metal which, in the UFF framework, is modelled with ε^UFF_Zn = 62.35 K compared to the Dreiding value of ε^Dre_Zn = 27.68 K. Then, depending on the type of framework atoms and their abundance within the unit cell, it should be possible to determine which MOFs will be sensitive to the choice of LJ FF before performing any molecular simulations. Indeed, Dokur and Keskin developed a simple factor, calculated using the difference in ε_i between LJ FFs and the number, n_i, of atoms i, to distinguish MOFs which are relatively insensitive, and MOFs which significantly depend on the choice of the LJ parameters (see supplementary note 6.2, ESI†).⁴⁵ Turning our attention now to some specific materials, we find that, similar to the previous scenario, often a particular feature of a MOF tends to magnify the expected deviations between LJ FFs. For example, in addition to OSOMIT, HELCUY, and GIMVAA, cluster I in Fig. 8(c) contains other Zn-metal MOFs such as BEPNIV (Fig. 9, top panel on the right) for which confinement effects cannot explain the observed sensitivities as these structures are quite open and feature relatively larger pores. These materials do however contain open metal sites (OMS).⁸⁵ As Zn is defined by a smaller LJ collision diameter in UFF, the closest distance of approach between adsorbing CO₂ and the OMS is shorter, which leads to an exponential increase in the electrostatic contributions to the interactions in these systems. This generally results in more sporadic adsorption behaviours using UFF parameters. From this observation we expect that large uncertainties are likely to be encountered when: (1) MOFs contain OMS that are accessible to the adsorbates; (2) the OMS belongs to a metal in which the difference in atomic diameters between LJ FFs is significant (e.g., Zn, Fe, Ti, Tc, and Ru); and (3) the partial charge assigned to the metal is sufficiently large. Caution is therefore warranted when interpreting the adsorption behaviours in MOFs with OMS using generic LJ FFs, particularly as they often fail to reproduce the adsorption data of their experimental counterparts.^85–89 Combined, the materials described here and in the previous scenario explain why poorer correlations between LJ FFs are observed using EQeq (Fig. 8(a and e)), and why for CO₂ in Fig. 8(b).

Finally, the third scenario is concerned with the accuracy of charge assignment using simplified schemes such as Qeq and EQeq. In many instances, charge equilibration methods do not capture the subtle differences in chemical bonding environments that are reflected by DDEC charges. Therefore, only modest agreement between DDEC and (E)Qeq charge schemes is observed over the entire CRAFTED-u database. A good example is given in Fig. S14 (ESI†), whereby (E)Qeq returns similar charges for all Zn atoms in all CRAFTED-u MOFs, while DDEC charges range between 0 and +1.5e (where e is the elementary charge). For specific materials, charge equilibration methods can fail spectacularly. In particular, we found that the Qeq scheme can return inappropriately high – and in some instances unphysical – charges. Such is the case for the MOFs populating the cluster of blue points (cluster II) below the parity line in Fig. 8(h). The spurious charges assigned in these materials arise for different reasons. For example, in EBOTOF (Fig. 9, lower panel on the right), Al metals are coordinated to highly electronegative atoms such as F, which represents a particular bonding environment for which Qeq fails. Similarly, abnormal Qeq charges are assigned in MOFs where Na, Ga, or In metals are present, such as in AYUWUM (Fig. 9, lower left panel). Ongari et al. noticed similar abnormalities in their study.⁴⁷ As to why Qeq fails in this regard, and EQeq does not, stems from the fact that Qeq uses a neutral charge-centre to express the Taylor series expansion in electric energy, while EQeq considers the Taylor series expansion around charge-positive metal cations. The neutral charge-centring is not able to accurately reproduce DDEC charges in alkali metals such as Na, but more generally it appears to fail for MOFs containing metals with a single electron in their outer orbital,⁴⁷ thus leading to similar deviations in In- and Ga-metal MOFs in the CRAFTED-u dataset. These outlier adsorption behaviours, which exist primarily for the Qeq scheme, are what contributes to the difference between Qeq and EQeq correlation lines in Fig. 8(b and f).

At this point, it is clear that the choice of the FF has a very strong impact on the adsorption behaviour, often leading to profound differences in the isotherms for the same material. These differences are likely to manifest in the material-level KPIs that depend on adsorption data from molecular simulations (see for example Fig. S16, ESI†), as well as in the performance estimates obtained using process simulations, discussed in the next section. Our results further emphasise that using simplified charge equilibration schemes in application to a heterogeneous database of materials is likely to result in unreliable predictions, which can further lead to significant scattering in the results and a lack of confidence in the rankings. From this analysis, we are now better equipped to establish a general connection between the different molecular FFs and the process performance distributions that arise from such modelling choices.

Quantifying process-level uncertainties emerging from the use of different molecular forcefields

In this section, we quantify the uncertainty between different FFs at the process-level. The objective here is two-fold: (1) we attempt to connect the material uncertainties described in the previous section to the process-level uncertainties in a systematic way; and (2) benchmarking protocols that quantify the impact of different FFs in conditions relevant to the desired application, e.g., carbon capture, will be required as improvements in molecular simulations continue to develop. Through our analysis, we aim to demonstrate the generality of our framework for evaluating future generations of molecular FFs.

For the above purposes, we make use of the metrics described in the Methodology section to map between distributions in adsorption isotherms and CO₂ purity-recovery Pareto fronts. A Pareto-dominance plot, constructed from for each molecular FF i, is displayed in Fig. 10(a). This plot demonstrates the tendency of a particular FF to dominate, or not dominate, relative to other FFs. In Fig. 10(b), the process-level uncertainties are visualised using the average hypervolume error, . Consistent with the discussion so far, we make comparisons amongst the different FFs by taking UFF + DDEC as the reference. For a complete comparison of the uncertainties between all pairwise FFs i and j, see supplementary note 7.1 (ESI†).

Fig. 10 Process-level uncertainty metrics for CO₂ purity-recovery optimisations. (a) Pareto-dominance plot constructed from values using the CRAFTED-u dataset. (b) Average hypervolume error, between i = UFF + DDEC and other molecular forcefields (j ≠ UFF + DDEC) for the CRAFTED-u dataset. Higher values of indicate a larger average uncertainty with the UFF + DDEC baseline.

Focusing on Fig. 10(a) first, it is clear that values constructed from the models based on EQeq charges are consistently lower than those using DDEC or Qeq charges. This means that the adsorption behaviours of materials generated using the EQeq charge scheme typically yield poorer process-level performances compared to DDEC and Qeq. Another interesting feature to note is the interaction between LJ FFs and the different charge schemes. Dreiding yields better average process performances than UFF when coupled with either DDEC or Qeq but poorer average process performances when coupled with EQeq. The molecular origin behind this inverse Pareto-dominance relationship is explained in supplementary note 7.2 (ESI†). Now, turning our attention to Fig. 10(b), we note a number of important trends. Relative to the UFF + DDEC baseline: (1) larger are typically observed when the charge scheme is changed rather than when the LJ FF is changed, i.e., ; (2) models based on EQeq charges are associated with larger compared to the models based on Qeq charges; and (3) the uncertainties associated with the choice of LJ FF is non-negligible. To expand on this final point, QOVWOO and XOVVIO are the 1st and 47th ranked materials according to the UFF + DDEC forcefield. Their performance curves deviate by a hypervolume of approximately 7%, which is equivalent to the average error between UFF and Dreiding in Fig. 10(b). The results presented here reasonably support the conclusions outlined in our discussion of the agreement between process-level rankings. In particular, they corroborate that the predictions are overall more sensitive to the choice of the charge scheme than to the choice of the LJ parameters, that the models based on EQeq charges agree less favourably with the DDEC baseline, and that the LJ FF can have a considerable impact on the process-level results.

On the contrary, we find that the connection between the process-level uncertainties and material-level correlations in Fig. 8(a, b, e and f) is not as straightforward. For example, given that a stronger correlation in the Henry's regime is observed between DDEC and EQeq in Fig. 8(b and f), one might expect that lower process-level uncertainties in Fig. 10(b) would be observed between DDEC and EQeq rather than DDEC and Qeq. However, different FFs often give rise to variations in several important features of both the CO₂ and N₂ adsorption isotherms simultaneously, such as the Henry's coefficient, local slope, nonlinearity, and total saturation capacity. The performance of an adsorbent is determined by an interplay of all of these features and by the competitive adsorption between CO₂ and N₂,^5,49–51,66 and so it is difficult to capture the process-level effects from a single metric such as the correlation between uptakes at low pressure. This makes establishing a general connection between molecular FFs and their process-level uncertainties a challenging task. It is possible to link individual materials to their process-level responses by a more in-depth analysis (for example, see supplementary note 7.2, ESI†), however the material-to-process mappings are in general nonlinear.

Exploring potential pathways towards uncertainty mitigation in multiscale HTCS workflows

In the previous sections we demonstrated that the typical choices and combinations of the FF parameters lead to inconsistent selection of the top performing materials and the relative rankings of the materials. We showed that both selection of the LJ parameters (UFF vs. Dreiding) and selection of the charge scheme has an impact on the predictions, however the inconsistency introduced by the charge scheme seems to be more significant compared to the choice of the LJ parameters. Furthermore, our results indicate that charge equilibration schemes may lead to unphysical results (Qeq) and strong deviations in the adsorption behaviours predicted by DDEC charges (both EQeq and Qeq).

These findings have significant implications for any further effort to implement fully in silico multiscale screening of porous materials. On the one hand, a strong statement would be that none of the combinations of the parameters are fully validated in experiments and, given the sensitivity of the process predictions to the selection of the FFs, it is unlikely that fully computational workflows will be able to produce reliable results. On the other hand, we argue that meaningful results can be obtained using this approach provided the parameter inputs and data outputs are handled with care. It is clear that shortcut methods such as (E)Qeq should not be used in HTCS studies, and so future screening studies should adopt ab initio levels of theory for partial charge assignment. The process performance of shortlisted materials identified using ab initio charges should therefore (ideally) be robust against perturbations in the LJ parameters. As a minimum, some estimate of the uncertainty introduced by different LJ parameter sets should be considered to build confidence in the predictions of these workflows. In either case, an accurate and comprehensive FF is the key bottleneck in computational screening, and we will return to this point in more detail in the next section. For now, we would like to explore two ideas to improve the consistency of predictions within the current state-ot-the-art.

The first idea is based on the hypothesis that amongst the CRAFTED MOFs, a group of materials exists whose adsorption behaviours are rather insensitive to the choice of the FF. It is further possible that among these materials, there are those that meet the required purity and recovery constraints and, in general, show good performance at the process level. We speculate that, given the insensitivity of the materials to the FF, it is likely that they will retain their performance in experiments as well. In this case, one can argue that the screening protocols could focus on searching for the materials in this class.

The second idea is as follows. We assert that ab initio charges are required for HTCS of MOFs for carbon capture. Unfortunately, it is still computationally challenging to employ these approaches for large databases of materials. Therefore, we would like to explore whether ML-based models^90–92 such as the message passing neural network (MPNN) model⁹⁰ and partial atomic charges in metal–organic frameworks (PACMOF) model⁹² can be used as a scalable but accurate alternative to the DDEC method.

Let us focus on these ideas in turn. To explore our initial hypothesis that there exists a subset of high performing MOFs that are relatively insensitive to the FF, we first define a taxonomy for describing CRAFTED MOFs by their ability to meet the CCS constraints for different FFs using four distinct consistency classes (Table 2). A small subset of 28 MOFs meets the CO₂ purity and recovery constraints on a consistent basis, regardless of the particular variation of the FF chosen. We denote these materials as being at the highest level of consistency: consistency level 1 (CL1).

Table 2 Summary of how CRAFTED MOFs are distributed amongst different performance consistency classes

Class^a	Label	Number of molecular forcefields (out of 6) for which a material meets the CCS constraints	Number of materials in class
a The lists of CSD codes for each of these classes are provided in supplementary note 8. b CL1 and CL2 class materials are differentiated from CL3 and CL4 on the basis of statistical significance. That is, the spread and skewness of the performance distributions for materials in these classes can be visualised using box-and-whisker plots which require a minimum of 5 measurements to construct.
CL1^b	Consistent CCS adsorbents	6	28
CL2^a	Semi-consistent CCS adsorbents	5	22
CL3	Inconsistent CCS adsorbents	[1–4]	168
CL4	Non-CCS adsorbents	0	472

---

The next question we pose is whether the high level of consistency in the CL1 group is because the adsorption behaviours in these materials are insensitive to the parameters of the FF or due to some other factors. To address this question, in Fig. 11(a) we visualise the relative ranking of CL1 class MOFs, using a box-and-whisker plot to represent the underlying distributions in process performances. It is evident from this subplot that the CL1 class contains materials which are characterised by both narrow distributions in performance, i.e., PEKKAS, and broad distributions in performance, i.e., NOQLOV. This suggests that a search for materials which consistently meet the CCS constraints is not sufficient to guarantee that their performance is insensitive to the choice of FF. Still, within the CL1 subset there are materials such as PEKKAS, QOVWOO, FIRVEH, and NAYZUK, whose narrow boxplots have the potential to confirm our hypothesis. By inspection of the individual adsorption behaviours and performance curves for select materials in Fig. 11(b), we observe a diversity of responses to the choice of FF. In every case however, some scattering in the adsorption data is observed. This means that the MOFs which are (rather fortuitously) characterised by narrow performance distributions achieve this level of consistency through different pathways. That is to say, the CO₂ purity-recovery Pareto fronts in these materials converge to similar high performing solutions (despite the scattering seen in the adsorption data) by using very different combinations of the PSA cycle parameters. Therefore, the results in Fig. 11 do not provide us with the type of evidence one would like to confirm our initial hypothesis. We arrive then to the conclusion that, in order to achieve consistent process performances in multiscale HTCS workflows for the right reasons, the inconsistencies which propagate from the material level must first be addressed.

Fig. 11 Relative ranking of CL1 adsorbents under uncertainty. (a) For each material, a box-and-whisker plot is used to visualise the spread and skewness of the process metrics used to evaluate CL1 class material performances. The process metric is the maximum CO₂ purity achievable (subject to the requirement that CO₂ recovery ≥0.9). The median process KPI across all forcefield combinations is used to rank materials, as indicated by the solid red line which is provided to guide the eye. Red crosses indicate outlier performance estimates. (b) Distribution of CO₂ purity-recovery optimisation results for select materials in class CL1. Top row shows the Pareto fronts for each material, middle row shows the CO₂ adsorption isotherms at 298 K, and bottom row shows the N₂ adsorption isotherms at 298 K. Solid lines and dashed lines indicate the UFF and Dreiding forcefields, respectively. Red, green, and blue colours indicate DDEC, EQeq, and Qeq charges, respectively. Note that the isotherms shown are generated from the dual-site Langmuir model which is used as input for the PSA process simulator.

This brings us onto our second idea. Here, we test whether the uncertainties which are introduced by charge equilibration schemes can be effectively mitigated by using the ML-based MPNN and PACMOF models. As a case study, we continue with the CL1 class of MOFs. These materials satisfy the CCS constraints across all FFs, and so we evaluate the accuracy of MPNN and PACMOF charges by optimising the energy penalty and productivity performances subject to the CCS constraints. This protocol is more robust than an unconstrained CO₂ purity-recovery optimisation as the energy-productivity process responses are more sensitive to changes in the adsorption behaviours.⁶⁶ The optimisations were first conducted using the ANN model and then refined using the PSA process model in order to calculate the hypervolume errors between Pareto fronts generated using different molecular FFs. The results are provided in Fig. 12, taking UFF + DDEC as the baseline for comparison.

Fig. 12 Process-level uncertainty metrics for energy penalty-productivity optimisations of CL1 class materials. Average hypervolume error, , between i = UFF + DDEC and other molecular forcefields (j ≠ UFF + DDEC). Higher values of indicate a larger average uncertainty with the UFF + DDEC baseline.

Disregarding the MPNN and PACMOF charges in Fig. 12 for a moment, we observe a similar response in the energy-productivity results to changes in the FF as the CO₂ purity-recovery results. However, it is clear that scattering in the adsorption data has a stronger influence on these KPIs. For instance, the average hypervolume error of 24% between LJ FFs (combined with DDEC charges) in Fig. 12 represents a 3-fold increase in the errors calculated at the CO₂ purity-recovery Pareto level (Fig. 10). While these results are determined for a reduced subset of MOFs, we generally expect an uncertainty of approximately 30% between LJ FFs using energy-productivity KPIs (see supplementary note 8.3, ESI†). Upon introducing the MPNN and PACMOF charges back into the discussion, the lowest uncertainties in Fig. 12 no longer occur between UFF and Dreiding using fixed DDEC charges. Instead, they occur between DDEC and MPNN/PACMOF charges using fixed UFF parameters, meaning that the prevailing sources of uncertainty identified previously (the choice of charge) has been effectively mitigated by modelling the electrostatics with ML-based charges. Interestingly, the uncertainty that is introduced by the MPNN and PACMOF models is comparable to the scattering seen between surrogate model predictions and PSA process simulations (Fig. 3). Efficient pre-screening of MOF databases could therefore be achieved by combining ML models for both the partial charge calculations and the PSA cycle optimisations without a significant loss of accuracy. Overall, our results suggest that MPNN and PACMOF can be used to assign charges in lieu of the DDEC method (see Fig. S18, S19 and S21, ESI†). This is further demonstrated in Fig. 13 whereby we see that, of the materials shown here, the adsorption behaviours and process-level predictions obtained using ML-based charges are in closer agreement with the DDEC baseline for all materials apart from NOQLOV.

Fig. 13 Distributions in energy-productivity Pareto fronts for different CL1 class materials. Top row shows the energy-productivity Pareto fronts for each material, middle row shows the CO₂ adsorption isotherms at 298 K, and bottom row shows the N₂ adsorption isotherms at 298 K. Note that the isotherms shown are generated from the dual-site Langmuir model which is used as input for the PSA process simulator. From left to right, each subplot shows a material with a lower degree of uncertainty between charge schemes.

Future prospects for HTCS of materials for carbon capture

While we strongly advocate for ab initio charges whenever possible, the results presented in Fig. 12 and 13 support our argument that ML-based charges are the next best option available to practitioners. Nonetheless, we recognise some of the shortcomings in our analysis (see supplementary note 9, ESI†) and understand that 28 MOFs cannot represent the diverse chemical landscape of all known MOFs. It is likely that expanding the scope of our study to more materials would expose some regions of poor agreement between ML-based and DDEC charges. Indeed, Liu and Luan³⁹ recently benchmarked several charge assignment methods and determined that charge predictions from PACMOF significantly deteriorated in MOFs which contained transition metals. However, as larger MOF databases with pre-computed ab initio charges begin to emerge,⁹³ we expect that improvements in the current ML models can be achieved by simply extending the training set to cover more materials with greater chemical diversity. Future efforts should therefore be directed towards developing more advanced MOF representations and ML architectures.

A problem still remains when we are posed with the question of which LJ FF to use. Both UFF and Dreiding are designed to be as generic as possible and thus are rather crude approximations of the complex interatomic interactions taking place during adsorption. These interactions are relatively weak compared to the energy variations in chemical bonds, and so the development of more reliable LJ FFs for adsorption in MOFs remains an open and important challenge. As with most areas of material science, ML has a possible role to play here. The interatomic potentials which describe short range interactions can be constructed using ML methods by approximating the high-dimensional potential energy surface (PES) using training data obtained from ab initio calculations.⁹⁴ The general and flexible nature of neural network models makes them a viable ML architecture for this purpose. Indeed, the so-called high-dimensional neural network potential (HDNNP) was introduced as early as 2007 by Behler and Parrinello,⁹⁵ and in principle considers the PES to be a sum of environment-dependent atomic energy contributions. HDNNPs can be computed many orders of magnitude faster than e.g., DFT calculations, they retain the accuracy of reference ab initio data, and can be scaled to large systems or time scales.⁹⁴

These concepts are only beginning to emerge in the field of MOFs. The first application of HDNNPs for MOFs was introduced by Eckhoff and Behler, where they predicted the negative thermal expansion of MOF-5 as well as its phonon density of states using HDNNPs.⁹⁶ Zheng et al. also recently developed a HDNNP to study to the chemisorption and diffusion of CO₂ in Mg-MOF-74.⁹⁷ While promising, this approach comes with its own set of open challenges. These include, but are not limited to, the construction of more informative feature vectors which can capture both the local chemical environments and MOF pore attributes, in developing an approach which can model material classes with a large number of chemical elements (the feature vector often scales with the number of chemical elements⁹⁸), and in generating the reference training data for MOFs with a large number of atoms in their unit cells.⁹⁹ Furthermore, extending this approach to predict adsorption in chemically diverse material databases would require transferable models: the training set must represent the problem at hand. A good example is the general-purpose ANI-2 potential,¹⁰⁰ which was trained on reference data for a large number of organic molecules and is thus capable of simulating systems containing (H, C, N, O, F, Cl, and S). One can envision a similar approach being adopted for mixed organic–inorganic complexes. An alternative, less data-intensive avenue would be to use active learning as a means to reliably explore the chemical space of MOFs.^99,101 In either case, a viable pathway towards developing a general-purpose, highly accurate ML potential for MOF screening applications exists provided that the appropriate training data, ML architecture, and feature vector can be combined. If such a method were to emerge, it has the potential to shift paradigms not only in the screening of MOFs for carbon capture, but in materials modelling in general.

Conclusions

In this multiscale computational screening study, we aimed to address the question: to what extent will the ranking of porous materials, and the selection of top performers identified in process modelling, depend on the choice of molecular forcefield? To do so, we generated distributions of CO₂ and N₂ adsorption isotherms in 690 MOFs using different forcefields (FFs) which represent the modelling choices commonly encountered in high-throughput computational screenings (HTCS) of materials for CO₂ capture. We then conducted a systematic uncertainty quantification study using PSA process-level modelling to determine one's ability to identify top-performing candidate materials consistently across different FF definitions.

Our results allow us to draw a number of general conclusions:

(i) Indeed, the computational ranking of materials appears to depend on the choice of the molecular FF to a significant extent.

(ii) From the pool of molecular modelling choices considered in this investigation, the choice of charge assignment scheme represents the largest source of uncertainty at the process-level.

(iii) Partial charges assigned by charge equilibration methods such as Qeq and EQeq are unable to reproduce the adsorption behaviours of charges derived from ab initio calculations such as DDEC with sufficient accuracy to guarantee consistent process-level rankings. As such, we recommend avoiding using charge equilibration methods when the electrostatic interactions are important for adsorption, such as CO₂ adsorption in MOFs. When ab initio charges are not available and are computationally unfeasible to obtain, ML-based charges are an attractive alternative that can effectively minimise the uncertainties originating from partial charge assignment. A simple extension in the pool of training materials can provide a quick remedy to the potential pitfalls of existing ML models, while future efforts should strive towards advanced MOF representations and machine learning architectures.

(iv) The process-level correlations between LJ FFs indicate that approximately 70% of the top performing MOFs may be mutually identified. This suggests that, in spite of considerable uncertainty, the search for MOFs using state-of-the-art multiscale HTCS is still more efficient than a random search. Nevertheless, we still lack a consistent, experimentally validated set of molecular parameters to describe the LJ interactions. Until this issue is addressed, one has to accept the considerable uncertainties embedded within process-level predictions which make use of adsorption data generated using generic LJ FFs.

We believe our work is an important step towards understanding the level of accuracy one can expect from multiscale screenings of materials for PSA-based carbon capture processes. The picture that emerges from our study is that, while HTCS remains a useful tool in the computational chemist's arsenal, it is not really possible to obtain material rankings with high resolution using this approach. In light of these observations, we see two ways to proceed with HTCS studies.

The first is to strive towards more consistent implementations of these workflows. As we mention above, the most direct route is to only use ab initio charges and, failing that, ML-based charges, to model the electrostatic interactions. Moving towards more consistent models of the short-range range dispersion/repulsion interactions would require two stages of implementation, we expect. In the short-term and until a truly universal methodology can be developed, opting for internally consistent LJ parameters within the community is an admirable pursuit. The rational choice would be to use the UFF forcefield, given its ability to simulate MOF databases with greater chemical diversities. However, we still recommend that the performance of shortlisted MOFs should be checked for robustness against perturbations in the LJ parameter sets. In the longer-term, we foresee great possibilities in the use of machine learned interatomic potentials as a means to simulate the short-range interactions in MOFs and, provided such a method can be developed, expect this approach to radicalise the way in which HTCS studies are conducted.

The second option is to alter the perceived utility of HTCS. Many technical barriers exist between identifying MOFs through HTCS and translating them into real, industrially viable separation materials. To name a few, a MOF must demonstrate good mechanical, thermal and moisture stability, low cost, and scalable synthesisability. A prime example is CALF-20¹⁰² which, despite its relatively modest process-level performance (see supplementary note 10, ESI†), is the only known MOF to satisfy enough of these technical requirements to succeed as a viable industrial-scale adsorbent. Therefore, rather than identifying the best performers via material rankings which we understand suffers from some reliability issues, one can instead try to identify good performers using HTCS and then shortlist materials based on all other factors which dictate the feasibility of an adsorbent. Given the importance of these material attributes, this approach could provide more meaning to the results obtained using HTCS and help expedite the transition between identifying promising MOFs and eventually tasking them to bench scale.

Author contributions

Conor Cleeton contributed to the data curation, conceptualisation, formal analysis, investigation, methodology, software, validation, visualisation, and writing – original draft. Felipe Lopes Oliveira contributed to the data curation, conceptualisation, methodology, and writing – review & editing. Rodrigo F. Neumann contributed to conceptualisation, methodology, and writing – review & editing. Amir H. Farmahini contributed to the conceptualisation, methodology, and writing – review and editing. Binquan Luan contributed to conceptualisation and writing – review & editing. Mathias Steiner contributed to conceptualisation and writing – review & editing. Lev Sarkisov contributed to the conceptualisation, formal analysis, methodology, validation, writing – original draft.

Code availability

All code created in this work is available on GitHub (https://github.com/ConorNCleeton/uncertainty_screening_mofs).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors would like to acknowledge the assistance given by Research IT and the use of the Computational Shared Facility at the University of Manchester.

References

1. A. Li, R. Bueno-Perez, D. Madden and D. Fairen-Jimenez, Chem. Sci., 2022, 13, 7990–8002.
2. A. H. Farmahini, S. Krishnamurthy, D. Friedrich, S. Brandani and L. Sarkisov, Chem. Rev., 2021, 121(17), 10666–10741.
3. S. Li, Y. G. Chung, C. M. Simon and R. Q. Snurr, J. Phys. Chem. Lett., 2017, 8, 6135–6141.
4. D. Wu, Q. Yang, C. Zhong, D. Liu, H. Huang, W. Zhang and G. Maurin, Langmuir, 2012, 28, 12094–12099.
5. A. K. Rajagopalan, A. M. Avila and A. Rajendran, Int. J. Greenhouse Gas Control, 2016, 46, 76–85.
6. T. D. Burns, K. N. Pai, S. G. Subraveti, S. P. Collins, M. Krykunov, A. Rajendran and T. K. Woo, Environ. Sci. Technol., 2020, 54, 4536–4544.
7. D. Yancy-Caballero, K. T. Leperi, B. J. Bucior, R. Richardson, T. Islamoglu, O. K. Farha, F. You and R. Q. Snurr, Mol. Syst. Des. Eng., 2020, 5, 1–26.
8. K. T. Leperi, D. Yancy-Caballero, R. Q. Snurr and F. You, Ind. Eng. Chem. Res., 2019, 58, 18241–18252.
9. A. H. Farmahini, D. Friedrich, S. Brandani and L. Sarkisov, Energy Environ. Sci., 2020, 13, 1018–1037.
10. S. Krishnamurthy, J. Boon, C. Grande, A. Lind, R. Blom, R. de Boer, H. Willemsen and G. de Scheemaker, Chem. Ing. Tech., 2021, 1–13.
11. S. G. Subraveti, K. N. Pai, A. K. Rajagopalan, N. S. Wilkins, A. Rajendran, A. Jayaraman and G. Alptekin, Appl. Energy, 2019, 254, 1–13.
12. D. Dubbeldam, K. S. Walton, T. J. H. Vlugt and S. Calero, Adv. Theory Simul., 2019, 2(11), 1900135.
13. F.-X. Coudert and A. H. Fuchs, Coord. Chem. Rev., 2016, 307, 211–236.
14. E. Haldoupis, J. Borycz, H. Shi, K. D. Vogiatzis, P. Bai, W. L. Queen, L. Gagliardi and J. I. Siepmann, J. Phys. Chem. C, 2015, 119, 16058–16071.
15. T. M. Becker, J. Heinen, D. Dubbeldam, L. C. Lin and T. J. H. Vlugt, J. Phys. Chem. C, 2017, 121, 4659–4673.
16. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909.
17. A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard and W. M. Skiff, J. Am. Chem. Soc., 1992, 10024–10035.
18. E. Braun, A. F. Zurhelle, W. Thijssen, S. K. Schnell, L. C. Lin, J. Kim, J. A. Thompson and B. Smit, Mol. Syst. Des. Eng., 2016, 1, 175–188.
19. L. C. Lin, A. H. Berger, R. L. Martin, J. Kim, J. A. Swisher, K. Jariwala, C. H. Rycroft, A. S. Bhown, M. W. Deem, M. Haranczyk and B. Smit, Nat. Mater., 2012, 11, 633–641.
20. N. S. Bobbitt and R. Q. Snurr, Mol. Simul., 2019, 45, 1069–1081.
21. M. J. Chiau Junior, Y. Wang, X. Wu and W. Cai, Int. J. Hydrogen Energy, 2020, 45, 27320–27330.
22. Y. J. Colón and R. Q. Snurr, Chem. Soc. Rev., 2014, 43, 5735–5749.
23. J. Stoll, J. Vrabec and H. Hasse, AIChE J., 2003, 49, 2187–2198.
24. S. Hamad, S. R. G. Balestra, R. Bueno-perez, A. R. Ruiz-salvador and M. Carlo, J. Solid State Chem., 2015, 223, 144–151.
25. C. I. Bayly, P. Cieplak, W. D. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
26. C. Breneman and K. Wilberg, J. Comput. Chem., 1989, 11, 361–373.
27. K. Sladekova, C. Campbell, C. Grant, A. J. Fletcher, J. R. B. Gomes and M. Jorge, Adsorption, 2020, 26, 663–685.
28. T. A. Manz and D. S. Sholl, J. Chem. Theor. Comput., 2010, 6, 2455–2468.
29. T. A. Manz and N. G. Limas, RSC Adv., 2016, 6, 47771–47801.
30. N. G. Limas and T. A. Manz, RSC Adv., 2016, 6, 45727–45747.
31. T. A. Manz, RSC Adv., 2017, 7, 45552–45581.
32. N. G. Limas and T. A. Manz, RSC Adv., 2018, 8, 2678–2707.
33. C. Campañá, B. Mussard and T. K. Woo, J. Chem. Theory Comput., 2009, 5, 2866–2878.
34. A. K. Rappe and W. A. Goddard, J. Phys. Chem., 1991, 95, 3358–3363.
35. C. E. Wilmer, K. C. Kim and R. Q. Snurr, J. Phys. Chem. Lett., 2012, 3, 2506–2511.
36. S. Ramachandran, T. G. Lenz, W. M. Skiff and A. K. Rappé, J. Phys. Chem., 1996, 100, 5898–5907.
37. B. A. Wells, C. De Bruin-Dickason and A. L. Chaffee, J. Phys. Chem. C, 2015, 119, 456–466.
38. M. Zhang and R. Fournier, J. Phys. Chem. A, 2009, 113, 3162–3170.
39. S. Liu and B. Luan, Nanoscale, 2022, 14, 9466–9473.
40. C. E. Wilmer, O. K. Farha, Y. S. Bae, J. T. Hupp and R. Q. Snurr, Energy Environ. Sci., 2012, 5, 9849–9856.
41. Z. Sumer and S. Keskin, Ind. Eng. Chem., 2016, 55, 10404–10419.
42. C. Altintas, G. Avci, H. Daglar, A. Nemati Vesali Azar, S. Velioglu, I. Erucar and S. Keskin, ACS Appl. Mater. Interfaces, 2018, 10, 17257–17268.
43. M. J. Lennox, M. Bound, A. Henley and E. Besley, Mol. Simul., 2017, 43, 828–837.
44. J. G. Mcdaniel, S. Li, E. Tylianakis, R. Q. Snurr and J. R. Schmidt, J. Phys. Chem. C, 2015, 119, 3143–3152.
45. H. Fang, H. Demir and D. S. Sholl, J. Mater. Chem. A, 2014, 2, 274–291.
46. C. Altintas and S. Keskin, Mol. Syst. Des. Eng., 2020, 5, 532–543.
47. D. Ongari, P. G. Boyd, O. Kadioglu, A. K. Mace, S. Keskin and B. Smit, J. Chem. Theory Comput., 2019, 15, 382–401.
48. D. Dokur and S. Keskin, Ind. Eng. Chem. Res., 2018, 57, 2298–2309.
49. A. K. Rajagopalan and A. Rajendran, Int. J. Greenhouse Gas Control, 2018, 78, 437–447.
50. K. T. Leperi, Y. G. Chung, F. You and R. Q. Snurr, ACS Sustainable Chem. Eng., 2019, 7, 11529–11539.
51. A. H. Farmahini, S. Krishnamurthy, D. Friedrich, S. Brandani and L. Sarkisov, Ind. Eng. Chem. Res., 2018, 57, 15491–15511.
52. Y. G. Chung, J. Camp, M. Haranczyk, B. J. Sikora, W. Bury, V. Krungleviciute, T. Yildirim, O. K. Farha, D. S. Sholl and R. Q. Snurr, Chem. Mater., 2014, 26, 6185–6192.
53. D. Nazarian, S. Camp and D. S. Sholl, Chem. Mater., 2016, 28, 785–793.
54. D. Dubbeldam, S. Calero, D. E. Ellis, R. Q. Snurr, D. Dubbeldam, S. Calero, D. E. Ellis and R. Q. Snurr, Mol. Simul., 2016, 42, 81–101.
55. danieleongari/EQeq: Charge equilibration method for crystal structures, https://github.com/danieleongari/EQeq, (accessed 21 July 2022).
56. F. L. Oliveira, C. Cleeton, R. Neumann Barros Ferreira, B. Luan, A. H. Farmahini, L. Sarkisov and M. Steiner, Sci. Data, 2023, 10, 1–12.
57. W. Lou, J. Li, W. Sun, Y. Hu, L. Wang, R. F. Neumann, M. Steiner, Z. Gu, B. Luan and Y. Zhang, Chem. Eng. J., 2023, 452, 139296.
58. Y. Hu, Y. Jiang, J. Li, L. Wang, M. Steiner, R. F. Neumann, B. Luan, Y. Zhang, Y. Hu, Y. Jiang, J. Li, L. Wang, Y. Zhang, M. Steiner, R. F. Neumann and B. Luan, Adv. Funct. Mater., 2023, 33, 2213915.
59. J. A. Ritter, K. C. Bumiller, K. J. Tynan and A. D. Ebner, Adsorption, 2019, 25, 1511–1523.
60. P. Nugent, Y. Belmabkhout, S. Burd, A. Cairns, R. Luebke, K. Forrest, T. Pham, S. Ma, B. Space, L. Wojtas, M. Eddaoudi and M. Zaworotko, Nature, 2013, 495, 80–84.
61. V. Benoit, R. S. Pillai, A. Orsi, P. Normand, H. Jobic, F. Nouar, P. Billemont, E. Bloch, S. Bourrelly, T. Devic, P. A. Wright, G. De Weireld, C. Serre and P. L. Llewellyn, J. Mater. Chem. A, 2016, 1383–1389.
62. Non-Linear Least-Squares Minimization and Curve-Fitting for Python—Non-Linear Least-Squares Minimization and Curve-Fitting for Python, https://lmfit.github.io/lmfit-py/, (accessed 5 December 2022).
63. A. L. Myers, AIChE J., 1983, 29, 691–693.
64. N. S. Wilkins and A. Rajendran, Adsorption, 2019, 25, 115–133.
65. S. J. Bhadra, A. D. Ebner and J. A. Ritter, Langmuir, 2012, 28, 6935–6941.
66. C. Cleeton, A. H. Farmahini and L. Sarkisov, Chem. Eng. J., 2022, 437, 135395.
67. M. Khurana and S. Farooq, Chem. Eng. Sci., 2016, 152, 507–515.
68. K. T. Leperi, R. Q. Snurr and F. You, Ind. Eng. Chem. Res., 2016, 55, 3338–3350.
69. J. C. Abanades, B. Arias, A. Lyngfelt, T. Mattisson, D. E. Wiley, H. Li, M. T. Ho, E. Mangano and S. Brandani, Int. J. Greenhouse Gas Control, 2015, 40, 126–166.
70. T. D. Burns, K. N. Pai, S. G. Subraveti, S. P. Collins, M. Krykunov, A. Rajendran and T. K. Woo, Environ. Sci. Technol., 2020, 54, 4536–4544.
71. K. Deb and H. Jain, IEEE Trans. Evol. Comput., 2014, 18, 577–601.
72. K. Deb, A. Pratap, S. Agarwal and T. Meyarivan, IEEE Trans. Evol. Comput., 2002, 6, 182–197.
73. L. Estupiñan Perez, P. Sarkar and A. Rajendran, Sep. Purif. Technol., 2019, 224, 553–563.
74. A. H. Farmahini, D. Friedrich, S. Brandani and L. Sarkisov, Energy Environ. Sci., 2020, 13, 1018–1037.
75. J. Beck, D. Friedrich, S. Brandani and E. S. Fraga, Comput. Chem. Eng., 2015, 82, 318–329.
76. M. T. M. Emmerich and A. H. Deutz, Nat. Comput., 2018, 17, 585–609.
77. Y. Cao, B. J. Smucker and T. J. Robinson, J. Stat. Plan. Inference, 2015, 160, 60–74.
78. X. Deng, W. Yang, S. Li, H. Liang, Z. Shi and Z. Qiao, Appl. Sci., 2020, 10, 569.
79. G. Avci, I. Erucar and S. Keskin, ACS Appl. Mater. Interfaces, 2020, 12, 41567–41579.
80. Y. Yan, Z. Shi, H. Li, L. Li, X. Yang, S. Li, H. Liang and Z. Qiao, Chem. Eng. J., 2022, 427, 131604.
81. S. Li, Y. G. Chung and R. Q. Snurr, Langmuir, 2016, 32, 10368–10376.
82. P. G. Boyd, A. Chidambaram, E. García-Díez, C. P. Ireland, T. D. Daff, R. Bounds, A. Gładysiak, P. Schouwink, S. M. Moosavi, M. M. Maroto-Valer, J. A. Reimer, J. A. R. Navarro, T. K. Woo, S. Garcia, K. C. Stylianou and B. Smit, Nature, 2019, 576, 253–256.
83. Y. G. Chung, D. A. Gómez-Gualdrón, P. Li, K. T. Leperi, P. Deria, H. Zhang, N. A. Vermeulen, J. F. Stoddart, F. You, J. T. Hupp, O. K. Farha and R. Q. Snurr, Sci. Adv., 2016 DOI:10.1126/sciadv.1600909.
84. J. A. Dunne, M. Rao, S. Sircar, R. J. Gorte and A. L. Myers, Langmuir, 1996, 12, 5896–5904.
85. M. Fischer, J. R. B. Gomes and M. Jorge, Mol. Simul., 2014, 40, 537–556.
86. G. Garberoglio, A. I. Skoulidas and J. K. Johnson, J. Phys. Chem. B, 2005, 109, 13094–13103.
87. S. Vandenbrande, T. Verstraelen, J. J. Gutiérrez-Sevillano, M. Waroquier and V. Van Speybroeck, J. Phys. Chem. C, 2017, 121, 25309–25322.
88. J. Zang, S. Nair and D. S. Sholl, J. Phys. Chem. C, 2013, 117, 7519–7525.
89. R. Mercado, B. Vlaisavljevich, L. C. Lin, K. Lee, Y. Lee, J. A. Mason, D. J. Xiao, M. I. Gonzalez, M. T. Kapelewski, J. B. Neaton and B. Smit, J. Phys. Chem. C, 2016, 120, 12590–12604.
90. A. Raza, A. Sturluson, C. M. Simon and X. Fern, J. Phys. Chem. C, 2020, 124, 19070–19082.
91. V. V. Korolev, A. Mitrofanov, E. I. Marchenko, N. N. Eremin, V. Tkachenko and S. N. Kalmykov, Chem. Mater., 2020, 32, 7822–7831.
92. S. Kancharlapalli, A. Gopalan, M. Haranczyk and R. Q. Snurr, J. Chem. Theory Comput., 2021, 17, 3052–3064.
93. J. Burner, J. Luo, A. White, A. Mirmiran, O. Kwon, P. G. Boyd, S. Maley, M. Gibaldi, S. Simrod, V. Ogden and T. K. Woo, Chem. Mater., 2023, 16, 54.
94. J. Behler, Chem. Rev., 2021, 121, 10037–10072.
95. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 146401.
96. M. Eckhoff and J. Behler, J. Chem. Theory Comput., 2019, 15, 3793–3809.
97. B. Zheng, F. L. Oliveira, R. Neumann, B. Ferreira, M. Steiner, H. Hamann, G. X. Gu and B. Luan, ACS Nano, 2023, 17, 5579–5587.
98. J. Behler, J. Chem. Phys., 2016, 145, 170901.
99. P. Friederich, F. Häse, J. Proppe and A. Aspuru-Guzik, Nat. Mater., 2021, 20, 750–761.
100. C. Devereux, J. S. Smith, K. K. Huddleston, K. Barros, R. Zubatyuk, O. Isayev and A. E. Roitberg, J. Chem. Theory Comput., 2020, 16, 37.
101. N. Bernstein, G. Csányi and V. L. Deringer, npj Comput. Mater., 2019, 5, 1–9.
102. J.-B. Lin, T. T. T. Nguyen, R. Vaidhyanathan, J. Burner, J. M. Taylor, H. Durekova, F. Akhtar, R. K. Mah, O. Ghaffari-Nik, S. Marx, N. Fylstra, S. S. Iremonger, K. W. Dawson, S. Partha, P. Hovington, A. Rajendran, T. K. Woo and G. K. H. Shimizua, Science, 2021, 374, 1464–1469.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ee00858d

---