Performance-Based Screening of Porous Materials for Carbon Capture

Computational screening methods have changed the way new materials and processes are discovered and designed. For adsorption-based gas separations and carbon capture, recent efforts have been directed toward the development of multiscale and performance-based screening workflows where we can go from the atomistic structure of an adsorbent to its equilibrium and transport properties at different scales, and eventually to its separation performance at the process level. The objective of this work is to review the current status of this new approach, discuss its potential and impact on the field of materials screening, and highlight the challenges that limit its application. We compile and introduce all the elements required for the development, implementation, and operation of multiscale workflows, hence providing a useful practical guide and a comprehensive source of reference to the scientific communities who work in this area. Our review includes information about available materials databases, state-of-the-art molecular simulation and process modeling tools, and a complete catalogue of data and parameters that are required at each stage of the multiscale screening. We thoroughly discuss the challenges associated with data availability, consistency of the models, and reproducibility of the data and, finally, propose new directions for the future of the field.


Introduction
Recent discoveries in material science and advances in computational chemistry are having a profound impact on the way we approach design and optimization of chemical processes, devices, and technologies. Traditionally, the workflow for the design of a process or a device would focus on a small number of materials available for experimentation and testing, as shown in the top panel of Figure 1. If performance of the material was not satisfactory, the experience gained in the process and the intuition of the investigator would guide the search for another material to be tried or suggest some modification of the existing material.
At the heart of an adsorption and membrane process is the material used as an adsorbent or a membrane. The efficiency of the process hinges on the characteristics of this material and the interplay between the material characteristics and process configuration. Recently, several new family of porous materials such as Metal-Organic Frameworks (MOFs), Zeolitic Imidazolate Frameworks (ZIFs), Covalent Organic Frameworks (COFs), Porous Organic Cages (POCs), Porous Aromatic Frameworks (PAFs), and polymers such as Porous Polymer Networks (PPNs) and Polymers with Intrinsic Microporosity (PIMs) have been discovered. A common motif, associated with these families, is a large number of (synthesized and hypothetical) members available within each family, tunability and exquisite control of structural characteristics of the materials such as surface area, pore size distribution (PSD) and surface chemistry. This prompted an extensive research effort to explore these new landscapes of structures to identify new porous materials with superior characteristics for adsorption applications, such as carbon capture The initial efforts in this field were led by the molecular simulation community, with various computational tools being used to obtain structural (e.g. surface area, porosity) and functional characteristics (e.g. equilibrium adsorption data) of these materials. These properties or metrics were then used to explore possible correlations between them and the function of the material in the actual application. An important question emerged from these early computational screening studies concerns the process descriptors or performance metrics: what descriptors and metrics should one actually adopt for ranking and selection of materials for a specific application? A useful metric must somehow reflect the essence of the process under consideration. For example, for methane storage the realistic metric is the working capacity, in other words the specific amount of methane released by the material when pressure is reduced from the storage pressure to the lowest pressure in the device, as oppose to the absolute capacity, corresponding to the lowest pressure being zero.
If for some applications, such as gas storage, a single metric may suffice the selection process, for other more complex, dynamic processes this is not possible. This was eloquently demonstrated by Rajagopalan et al. [6] by comparing a broad range of traditional and new separation performance metrics developed over the years and the actual performance of the material in the process simulator using post-combustion CO2 capture as a case study.
In fact, a significant amount of literature and studies have been accumulated over the years on design and optimization of pressure, vacuum, temperature, concentration, electric, and microwave swing adsorption processes, from simplified equilibrium models to more advanced numerical approaches [7][8][9][10][11][12][13][14][15][16][17][18]. Typically, these studies focus on a particular process configuration and conditions, while the cycle configuration is optimized to meet specific process targets. In the case of the aforementioned carbon capture application, the targets (or constrains of the process) are 90% recovery of the CO2 from the feed with 95% purity, as recommended by the DOE based on the emission control targets and storage requirements [19]. The efficiency of the process and hence performance of the material for the process is then assessed from the perspective of two metrics: productivity, in other words the amount of CO2 captured per unit of time by a unit of volume of the adsorbent, and energy penalty, which is the energy required to capture a mole of CO2 in the process. These two metrics are in competition with each other and a complex trade-off between them cannot be captured using simplified equilibrium-based figures of merits.
The co-current developments in computational screenings based on molecular simulations and in advanced process simulations invariably led to the following proposition: what if the screening of porous materials for dynamic adsorption processes can be implemented using realistic process simulations while the microscale properties of materials are provided by molecular simulations? This multiscale screening protocol is schematically depicted in Figure 2. According to this diagram, molecular simulations can be used to obtain equilibrium data (e.g. adsorption isotherms), dynamic properties (e.g. micropore diffusivity) or other materials characteristics (e.g. thermal properties), if needed. This information is then fed into a process simulator and the performance of the materials is assessed using the metrics previously developed for dynamic adsorption process analysis. The first examples of such a multiscale approach were published in two pioneering studies by Hasan et al. [20] for in silico screening of zeolite materials in the context of carbon capture, and by Banu et al. [21] for computational screening of MOFs for hydrogen purification. The early endeavours into the field of multiscale approaches to material screening also exposed a number of challenges. These challenges are associated with consistent and reliable transfer of data and information between the different levels of the simulation (from molecular simulations to process simulations), sensitivity of the process simulation predictions to the properties that cannot be obtained from molecular simulations, lack of experimental validation of the process simulation predictions, the accuracy of the produced material rankings, and lack of accurate molecular force fields for new families of materials, just to name a few.
It is also clear that multiscale approaches where one is able to seamlessly progress from a material structure (virtual or real) to the performance of the material in the actual process or device will become of immense importance in the nearest future. With the advent of machine learning and Quantum-Mechanical methods we are witnessing the dawn of material-driven process design, which will have a profound impact on a number of technologies and applications.
Hence, this review is prompted by recognition of the importance of this emerging field for multiscale material screening and discovery, and challenges that have been already encountered in the early studies. Here, we seek to provide a single source of information for scientists working across various fields for the development of multiscale materials screening approaches for adsorption-based postcombustion carbon capture. This includes chemists and materials scientists working on the development and characterization of new adsorbents; computational chemists and molecular modellers who develop atomic force fields and new computational tools for molecular simulations; and finally experts in the field of process modelling and optimization. We note that although this review deliberately focuses on the post-combustion carbon capture using PSA/VSA processes, the multi-scale workflow developed for this purpose and the challenges associated with advancement of this approach will be similar for a wide range of other separations processes such as hydrogen separation, oxygen purification, air separation and so on.
Throughout this review, we aim to highlight the fact that development of accurate and efficient multiscale workflows for realistic screening of porous materials can only be successful, if scientists working on different elements of these workflows are aware of the requirements of other parts. We also hope that the current review can encourage cross-disciplinary collaborations in this emerging field and lead to the development of multiscale screening tools to be used in a variety of settings, from chemistry labs to chemical engineering pilot plants. With this in mind, the specific objectives of this article are as follows: i. Review recent major developments in the field of multiscale approaches to material screening for adsorption separation processes with a focus on post-combustion carbon capture.
ii. Provide a single source of information on the basic principles of molecular and process simulations, data required at each stage, sources of data, and sources of uncertainties.
iii. Review the key challenges in the implementation of the multiscale screening strategies and how they can be tackled. iv. Outline future developments and trends in this emerging field.
To a significant extent this review is motivated by many conversations we had with our chemist and material scientists' colleagues who, having synthesized or proposed a new material, wished to know what it would do in the actual application. Testing it in a multiscale workflow would be the answer to this question. With the audience specified above, our philosophy has been to introduce the required elements of the multiscale process workflows at the level accessible to non-specialists, while the process and molecular simulation communities have more in-depth reviews and sources available to them.
To this end, the review is divided into eight main sections. After this introduction, Sections 2, 3, 4 and 5 will cover the application in question (post-combustion carbon capture), explain basic elements of pressure/vacuum swing adsorption processes, discuss hierarchy of metrics that can be used for selection and screening of porous materials for this application, and provide a historical perspective on how computational screening methods evolved over the last 10 years towards current multiscale workflows. These four sections are intended for readers who want to briefly explore recent developments in this field and are familiar with the key drivers, contributors and contributions in the field. Section 6 mirrors in its structure the multiscale workflow depicted in Figure 2. Here, we will cover practical aspects associated with material databases and the tools available for structural characterization of materials that are currently collected by these databases. Next, we will move to introduce the fundamentals of molecular simulations and process modelling. We will explain how these elements should be used together and as part of a multiscale workflow for materials screening. For each method, we will also introduce available simulation tools and software packages that can be used for performing these types of simulations. Our emphasis will be on explaining what data are required at each stage and what information is obtained at each level, but we will also discuss the gaps in the methods that need to be addressed. In Section 7, we will explore current challenges in the field of multiscale materials screening and will provide our suggestions for addressing them, which we hope will stimulate further cross-disciplinary approaches and collaborations. Finally, in Section 8, we finish the review with a brief discussion on future opportunities and possible directions of the field. . Adapted from Reference [4].
Carbon capture and sequestration (CCS) remains one of the key priorities in addressing the global climate change. This is the area where additional energy penalty associated with preventing carbon dioxide emission from power plants is the most significant barrier to the implementation of CCS technology, and any advance in this domain will likely have a profound impact on our ability to control carbon dioxide levels. For this reason, CCS has been one of the most explored applications in the context of computational screening of new materials, such as zeolites, MOFs, ZIFs and others. This is also the area where the multiscale workflow approaches have made the most significant progress. Hence, CCS and in particular post-combustion capture is a logical focus of this review.
Given the intended target audience of this review, it is useful to introduce the basic concepts of postcombustion carbon capture, while referring the interested reader to the more specialized and extensive sources on the topic [22][23][24][25][26][27][28].
The 2005 IPCC [4] committee identified three possible technologies for carbon capture from power plants, the most significant stationary CO2 emitter globally: pre-combustion carbon capture, oxy-fuel process and post-combustion carbon capture (Figure 3 a). In pre-combustion capture fuel reacts with oxygen (or air) and steam. This produces so-called syngas (synthesis gas) composed predominantly of carbon monoxide and hydrogen. In the water-shift reactor, this mixture reacts with steam to produce carbon dioxide and more hydrogen. Carbon dioxide is then separated from the mixture and the remaining purified hydrogen is used as a clean fuel in various processes. The idea of the oxyfuel process is to use pure oxygen for combustion. This oxygen is produced in the air separation step, which naturally comes with energy cost. However, as the process produces pure carbon dioxide during the combustion step, it does not require any carbon dioxide separation step, saving the costs down the line. Finally, in the post-combustion process carbon dioxide separation is applied to the flue gas from a standard power plant.
Post-combustion capture is the only technology that can be retrofitted onto the existing power plants and therefore is a promising approach in short and medium terms. In fact, detailed analysis of the US National Energy Technology Laboratory's (NETL) CCS database shows that there are currently more than 30 active post-combustion carbon capture plants around the world [29]. This is illustrated in Figure 4. In addition, post-combustion capture can be applied to hard to decarbonise emissions such as from industrial processes and to power plants converted to bioenergy (BECCS) which would enable negative emissions.  [29].
The composition of the flue gas is typically 15-16 vol% CO2, 5-7 vol% H2O, 3-4 vol% O2, and 70-75 vol% N2 for coal-fired power plants. In addition, the flue gases may contain trace amounts (tens and hundreds of ppm) of carbon monoxide, SOx and NOx. This stream is at 1 bar and 50-75°C [30]. We note, however, that most of the design efforts focus on a simplified separation operation involving only a binary mixture of CO2 and N2 at 1 bar and temperatures below 40°C.
A viable carbon capture technology must remove 90% of carbon dioxide from this flue gas and produce it with 95% purity [19]. The purity constraint is mostly dictated by the requirement to compress the product CO2 gas to 150 bar for further transportation and geological storage. Higher proportion of nitrogen in this stream would incur higher compression costs. These targets set the basis for the comparison of technologies proposed for this task.
Traditional approaches for carbon capture from power plant streams employ solvent-based (e.g. amine) absorption processes. It is estimated that the best absorption technologies incur a parasitic energy penalty of about 1.3 MJ per kg CO2 captured [31]. This is associated with a significant energy demand of the solvent regeneration step. Any new technology proposed for carbon capture must demonstrate that it is economically more viable (i.e. has lower energy penalty) than the reference, state-of-the-art amine absorption processes.

Pressure and Vacuum Swing Adsorption for Post-Combustion
Carbon Capture Figure 5. A schematic 4-step VSA cycle for separation of CO2 and N2 (a), difference of PSA/VSA/TSA processes illustrated using equilibrium adsorption isotherms of CO2 and N2 (b).
The main objective of this section is to introduce the key concepts and terminology associated with the pressure/vacuum and temperature swing adsorption processes that are required later in the article. The essential principle behind adsorption separation is that the components of the gas or liquid mixture somehow interact differently with the porous material and this difference can be exploited to separate them. Depending on the nature of this difference, we can distinguish i) kinetic separations, in which diffusion of molecules of the mixture in and out of the material happens at a significantly different rate; ii) molecular sieving, where one of the components of the mixture is simply too bulky to fit in the pores of the structure; or iii) equilibrium separations, where one of the components interacts more strongly with the porous structure via intermolecular interactions, which is the largest group if industrial processes. The PSA/VSA processes under consideration in this article belong to this last equilibrium group of separation processes.
To illustrate the principles of a PSA/VSA processes, let us consider the diagram depicted in Figure 5 (a), which shows different phases of a typical PSA/VSA process. The main element of this diagram is the adsorption column (schematically shown as just a rectangular box) filled with the porous material, or adsorbent. In the first step (adsorption) the feed is introduced in the column. Stronger interacting components (called heavy components) are preferentially adsorbed by the porous material in the column, changing the composition of the gas phase. As a result, the product gas stream leaving the column on the other side (so called, raffinate) is rich in the light components (weakly adsorbing components of the mixture). At some point in time, the column becomes saturated and will not be able to adsorb anymore of the heavy components. At this point, the adsorption step should be stopped, and the column should go through the regeneration or desorption phase. This phase may consists of a preliminary pressure reduction step (the blowdown step) followed by further reduction of pressure (the evacuation or extraction step), moving the process to the conditions associated with the low loadings on the isotherm and causing desorption of the heavy component (Figure 5 (b)). The column is then re-pressurized and goes through the adsorption step again.
As can be seen from the simplistic description above, the PSA/VSA process is a cyclic process, where the basic unit of the process, the adsorption column, goes through the repeating phases of adsorption and desorption. In the example above, we used pressure change (or swing) on the adsorption isotherm to regenerate the column as depicted in Figure 5 (b). Alternatively, we could have used higher temperature for regeneration. Indeed, as adsorption from the gas phase is an exothermic process, a higher temperature will shift the equilibrium to lower loadings, leading to desorption. This process is called Temperature Swing Adsorption (TSA). A combination of Pressure and Temperature Swing is also possible (PTSA) and the trajectory of conditions associated with this process is also shown in Figure 5 (b).
The difference in the equilibrium amount adsorbed between the adsorption and desorption cycle is called the working capacity. If the PSA system is cycling between ambient pressure and vacuum, then it is called a Vacuum Swing Adsorption (VSA) process. The main additional energy cost of PSA/VSA processes is associated with pulling the vacuum (VSA) and compression (PSA). Hence, the work of vacuum pumps and compressors become a key ingredient in the assessment of economic viability of the PSA/VSA processes.
For the PSA/VSA adsorption process to operate continuously, the actual plant consists of several columns going through various stages of the cycle. The number of units and how they are arranged is called process configuration. The types of steps involved, the timing of the steps within a single cycle, their duration and other parameters constitute a cycle configuration. Developing process and cycle configurations, in order to lower energy penalty and increase productivity constitute the main objective of the PSA/VSA design process.
In case of the post-combustion separation process, carbon dioxide is the heavy component and nitrogen is the light component. Zeolite 13X is the most explored material for this application, both in process modelling and in pilot plant studies. This material is hydroscopic and will adsorb water present in the flue gas, leading to higher cost of the process. We will discuss the current state-of-the-art technology and the associated energy costs in various process configurations in Section 6.3. The brief introduction provided in this section serves only to establish the most essential elements of the PSA/VSA processes, while for the more extensive reviews of gas adsorption separation processes for carbon capture the reader is referred to more specialized and extensive sources [27,[32][33][34][35].

Hierarchy of Performance Metrics for Materials Screening
In Section 2 we described the problem in hand: to capture CO2 from flue gas of a power plant with 90% recovery and 95% purity. Imagine now that we intend to identify the best adsorbent material for this task using a PSA/VSA process (as briefly described in Section 3) from a cloud of many thousands of possible porous materials. To do so, we need a suitable performance indicator (i.e. metric) which can correctly quantify separation performance of porous materials and sufficiently discriminate between performances of similar materials. Hitherto, a large number of performance indicators has been proposed for the assessment of materials separation performance. In this section, we review the most important of these indicators as reported in the literature, while focusing predominantly on their nature, classification, and availability. The information provided here will form the basis of the following discussion in the next section where we will illustrate how application of these metrics in the field of computational materials screenings evolved over the years leading to wider adaptation of process-level metrics for materials ranking.
Colloquially speaking, one would want to select the best material for a particular application simply by looking at its structure. In more scientific terms, the first group of metrics contains descriptors of the material structure: its porosity, density, surface area, pore size distribution (PSD) and so on [36,37]. These properties can be either obtained from the experiments, as a part of the standard characterization procedure for every new material synthesized, or from the computational methods as will be discussed later in this article. We call this group of metrics intrinsic structural material metrics (ISMM). These structural metrics do not tell us anything about how material interacts with its environment. Functional behaviour of the material is described by adsorption equilibrium data (adsorption isotherms, Henry's constants of adsorption, adsorption capacity), transport characteristics (diffusivity), thermal properties, etc [38][39][40][41]. These properties can be termed intrinsic functional material metrics (IFMM).
In separation applications, adsorption is a competitive process between two or more adsorbing species. Naturally, to characterize this competition we need a metric which would compare the behaviour of the material with respect to competing species. For example, selectivity is the ratio of loadings for two gases and at low pressure can be expressed simply as the ratio of the two Henry's constants. Selectivity is the simplest metric from the group of hybrid metrics (HMM), which combine various adsorbent metrics mentioned above, to more accurately discriminate between adsorbents with different separation performances. Examples of these metrics include: adsorption figure of merits (AFM) [42], sorbent selection parameter (SSP) [43], separation factor (SF) [44], adsorbent performance indicator (API) [45], and adsorbent performance score (APS) [46]. Mathematical definitions of these metrics are provided in Table 1.
One important step in the progress in the development of more realistic metrics for materials screening was the realization that selectivity and working capacity are not necessarily representative of the economic drivers of gas separation processes. To address this limitation, new screening metrics were developed to exploit the correlations between adsorption characteristics of porous materials and the economic drivers of separation processes. The first prominent example of such evaluation metrics is the separation performance metric (SPP) developed by Braun et al. [47]. SPP was developed to be representative of the most important economic drivers for separation of CO2 from natural gas mixtures [47]. It assumes equilibrium adsorption and desorption in the PSA/TSA/PTSA processes in order to calculate the value of an objective function, which accounts for the amount of captured target gas (e.g. CH4), amount of adsorbent material used, and the total energy required for the separation process [47]. The assumption of a process performing fully under equilibrium represents an ideal case scenario, however this condition is not always achieved in dynamic separation processes such as PSA/VSA. The other limitation of SPP metric is that instead of using conventional cost indicators (e.g. capital and operating costs), SPP assumes that all process costs scale with the amount of adsorbent ( ) used in the separation unit [47]. However, there are cases where a significant portion of the capital costs is independent of the amount of material used in the process. If these contributions of the capital cost become significantly larger, the amount of material used in the separation unit will become irrelevant [47]. Comparison of SPP, SSP and API metrics with detailed process modelling indicates that for CO2/CH4 separation, SPP surpasses the other two evaluation metrics in terms of accuracy [47].
Another important example of new evaluation metrics is the parasitic energy (PE) which was first used by Lin et al. [48] and Huck et al. [49] for evaluation of different classes or porous materials for postcombustion carbon capture. In their analysis, the additional energy required for adsorption carbon capture process consists of: (1) energy to heat the adsorbent material, (2) energy to supply the heat of desorption which is equal to the heat of adsorption, and (3) energy needed to compress CO2 to 150 bar which is a standard requirement for transport and storage [48]. Based on this, the authors formulated a simplified expression for the parasitic energy of a CCS process as a combination of the thermal energy requirement and the compressor work [48]. In the definition of parasitic energy provided by Lin et al. [48] equilibrium adsorption and desorption is assumed. As mentioned before, this may not be always the case in dynamic PSA/VSA systems. The parasitic energy curve obtained for screening of a large group of porous materials is shown to be a good reference to benchmark other newly discovered materials [48].
Inadequacy of screening metrics which are solely linked to the adsorbent properties and not their performance at the process level has been recently demonstrated by Rajagopalan et al. [6] using a case study for post-combustion CO2 capture. Without intending to repeat the entire argument here, one may consider an example selectivity of a candidate material for CO2/N2 separation using PSA process. On its own, a high value of selectivity is unlikely to be enough to select the material for CO2 separation. For instance, if the material has very low capacity the operation is likely to be very costly, despite high selectivity of the material. This study clearly demonstrates that for complex, dynamic adsorption processes such as PSA/VSA processes for carbon capture, the realistic performance of a specific material must be assessed in the actual process, by performing process simulation and optimization under realistic conditions. For this purpose, a new class of evaluation metrics is required. The metrics used to assess performance of porous materials at the process level are therefore called process-level metrics (PLM) in this review. In this case, a trade-off curve between overall energy penalty of the process and its productivity is used as an evaluation metric for materials screening [6,50,51]. This metric is normally calculated under certain constrains for purity and recovery of the product in the PSA/VSA processes. For example, for the case of post-combustion carbon capture, the minimum purity and recovery of CO2 are maintained at 95% and 90%, respectively. Energy penalty and productivity not only are more realistic measures of process performance (because they are obtained from detailed process modelling and optimization), they are also more directly related to the economic drivers of the separation process. Therefore, the next natural step in developing realistic evaluation metrics for materials screening is to link the existing process modelling platforms to technoeconomic analyses of the process because the ultimate goal of any separation unit is to achieve the design objective at the lowest cost [52][53][54]. Khurana and Farooq have extended this concept to include a comprehensive costing framework for the entire carbon capture plant [52,53]. Their integrated optimization framework looks at the separation cost in terms of $/tonne of CO2 captured or $/tonne of CO2 avoided, where the latter is defined as the difference between emissions of two power plants, one without a capture unit and the other with a capture unit but both producing the same net amount of electricity [52]. Fully integrated techno-economic analysis of carbon capture plants or any other industrial separation facility can be a daunting task for the purpose of screening of large groups of adsorbent materials that are currently available. As a result of this limitation, more recent studies have attempted to develop general evaluation metrics (GEM) that are strongly correlated with the results of the detailed techno-economic analyses [55]. Usually, this is achieved by combining all previously known evaluation metrics into a more general one (i.e. GEM) and then reducing complexity of the GEM by removing the elements whose contribution into the correlation coefficient is insignificant [55,56]. Leperi et al. [55] has shown that this approach is quite promising for the development of universal screening metrics that simultaneously take into account most important characteristics of the process associated with adsorbent material, process optimization and overall economic cost of the plant. Development of new GEMs can particularly benefit from recent advances in machine-learning techniques, if adequately large datasets of techno-economic materials were available for training the GEM function.
From the provided review of the hierarchy of metrics, one could make an impression that if the most accurate assessment of the material performance is provided by detailed process and plant models, then this should be the standard level of description in all material screening protocols. This, however, does not take into account, the computational cost of these metrics. Once the equilibrium adsorption data are available, the hybrid methods provide effectively an instant assessment of the material candidate. Process simulation of a single cycle configuration for a PSA/VSA process may be done in a few minutes in conventional 4-cores CPUs, whereas cycle optimization for the best performance may take tens of hours. This computational price tag applied to thousands and tens of thousands of materials would still make routine use of screening of all materials at the process level unaffordable. Hence, this is still an ongoing area of research to develop a multistage screening process, where efficiency of process optimization are improved using novel numerical techniques, or alternatively some preliminary screening is done using hybrid metrics/simplified process models, while accurate process modelling and optimization is only carried out for a selected group of promising materials.   [47]. For PE, Q, η and Wcomp are the thermal energy requirement, Carnot efficiency and compressor work respectively. For GEM, WCmod stands for modified working capacity as defined in the corresponding reference [55].

Computational Screening of Porous Materials: A Historical Perspective
In the previous section, we discussed what metrics are available for material screening in adsorption applications through the prism of metric hierarchy from very simple "intrinsic" metrics to processlevel metrics. In this section, we take a different, historical perspective on the development of computational screening strategies. This perspective will allow us to review how this field has evolved over time towards current multiscale workflows that incorporate elements of different type simulation techniques and performance indicators.
The first material screenings can be tracked back to more than 10 years ago where this terminology started to be used more widely by the scientific community [58][59][60]. In a pioneer study published in 2010 [38], Krishna and van Baten employed the configurational-biased Monte Carlo (CBMC) and molecular dynamics (MD) simulations to examine adsorption, diffusion, and permeation selectivities for separation of CO2/H2, CO2/CH4, CO2/N2, CH4/N2 and CH4/H2 mixtures in a number of zeolite, MOF, ZIF and carbon nanotube (CNT) structures. Their studies provided useful guidelines to the optimum choice of microporous layers that should be used in membrane separations representing a compromise between selectivity permeability and the permeability itself. This study also emphasized the importance of correlations between pore space properties (pore volume, limiting pore diameter, etc.) and transport properties (e.g. diffusion and permeation) in these classes of porous materials.
Building on the importance of the pore structure characterization, Haldoupis et al. [37] analysed pore sizes of more than 250,000 hypothetical silica zeolites to compute the size of the largest adsorbing cavity and pore-limiting diameter for all zeolites. This information can be used to reveal the range of adsorbate molecules that can possibly diffuse through each zeolite. Additionally, the authors computed Henry's constant of adsorption and diffusion activation energy for CH4 and H2 for a subset of 8000 zeolites. From the diffusion activation energies, they were able to estimate diffusivity of each adsorbate using the Transition State Theory (TST). In this study, they demonstrated that using a combination of molecular simulation techniques, one can rapidly assess adsorption properties of a large group of nanoporous crystalline materials for a particular application [37].
Application of computational materials screening approaches took another step forward in 2012 when two major studies were published. Namely, Snurr and co-workers used a library of 102 building blocks and a "tinker-toy" algorithm to assemble a database of 137,953 hypothetical MOFs [61]. Using geometric characterization tools and Monte Carlo simulations, they explored their database to identify the most promising structures for methane storage. From this perspective, this is the first example of a computational screening strategy applied to a large group of MOF materials. Later in the same year, Snurr and co-workers [62] simulated adsorption of CO2, CH4, and N2 in more than 130,000 hypothetical MOFs from the same database and subsequently examined their potential for CO2 capture using five different performance indicators (i.e. metrics) including CO2 uptake, working capacity, regenerability, adsorption selectivity and sorbent selection parameter (as defined in Table  1). They showed that although the resulting structure-property relationships between pore size, surface area, pore volume, and chemical functionality provide several leads for design of new porous materials, none of the above metrics is actually a perfect predictor of CO2 separation performance. The studies of Snurr and co-workers introduced several concepts that are now central to the computational screening strategies of porous materials. The concepts can be formulated as follows: i. The modular nature of MOFs allows the use of simple, tinker-toy algorithms to assemble new hypothetical structures simply by linking the building blocks along the appropriate topology. This idea can be extended to other new classes of materials (ZIFs, COFs, etc).
ii. Each material within the database can be explored in terms of structural properties (surface area, pore volume, largest pore size) and functional properties (capacity with respect to a particular gas under specific conditions, Henry's law constant of adsorption, selectivity). These properties can be used to classify, compare, and organize materials within the database.
iii. Computational screening studies calculate properties mentioned above. Two or more properties correlated to each other form clouds of data points, which can be explored to reveal some promising structure-property relations.
As will be discussed in the following, further studies in this emerging field also identified some deficiencies of the original database of hypothetical MOFs constructed by Wilmer et al. [61] and, consequently, identified challenges and new directions of research. These can be summarized as follows: i. Structures assembled in the tinker-toy algorithms require further accurate structure optimization using Quantum Mechanical (QM) methods to be more realistic.
ii. In general, we need systematic approaches to organize structures into databases that can be used in molecular simulations.
iii. Early computational and experimental studies in 2000-2010 very clearly demonstrated that the  accurate molecular force fields are lacking for new classes of porous materials interacting with gases  and liquids. A particularly striking manifestation of this was the failure of the conventional force fields  to describe interaction of MOFs featuring open metal sites with carbon dioxide or unsaturated hydrocarbons. Interaction of adsorbents with water presents wholly separate and substantial challenge. This prompted the simulation community to put significant effort into the development of a new generation of force fields based on accurate QM potential energy surface. Despite some initial and significant advances, atomic force fields remain mainly specialized and non-transferable; and this is still very much a remaining challenge and an ongoing area of research.
iv. Early studies would use several simple, well-known algorithms to obtain structural characteristics of the porous material. Later a number of comprehensive and versatile tools were developed (Zeo++ [63], Poreblazer [64], ZEOMICS/MOFomics [65,66]) to calculate geometric descriptors of porous materials. These tools vary in the algorithms, characteristics they calculate and in the philosophy of use and access.
v. Development of Machine Learning algorithms is needed to establish structure-property relationship within the databases and drive the discovery of new materials with desired functionalities.
At the same time, Smit and co-workers [48] also published a new study on screening of hundreds of thousands of zeolite and ZIF structures using the parasitic energy (PE) as a promising metric for evaluation of materials performance in the context of post-combustion carbon capture. PE is defined as the minimum electric load imposed on a power-plant by a carbon capture and storage (CCS) unit. At the molecular simulation level, a combination of grand canonical Monte Carlo (GCMC) simulation, energy grid construction method and Widom test particle insertion technique were employed to obtain equilibrium adsorption characteristics of materials. The PE metric was then used to search for materials that have the potential to reduce the parasitic energy by 30-40% compared to the conventional amine-based absorption technologies [48]. This study proposed a theoretical limit for the minimal parasitic energy that can be achieved for a particular class of porous materials.
A series of early articles by Sholl and co-workers [60,67], and Keskin and colleagues [68][69][70] had laid the foundation of the computational screening methods for membrane gas separations. This was followed by Kim et al. [71] publishing a major study on screening of over 87,000 different zeolite structures for permeation separations [71]. In this publication, the authors estimated the diffusion coefficients of CO2, N2 and CH4 using free energy calculations and the Transition-State Theory (TST) and identified general characteristics of the best-performing structures for CO2/CH4 and CO2/N2 membrane separations. For CO2/CH4 separation, they predicted a structure that outperformed the best known zeolite structure by a factor of 4-7 based on the required area of an ideal membrane which is shown to be mainly dominated by and inversely proportional to the CO2 permeability in the system [71]. In comparison with the results of Haldoupis et al. [37], Kim et al. demonstrated that screening of porous materials based on purely geometric approaches may deviate from what is predicted from a more advanced energy-based analysis [71].
The above study was followed by two other publications from the same group with a greater emphasis on MOFs as an emerging group of porous solids for adsorption separation applications. The first study was published in 2014 by Sun et al. [39] where 12 materials including six MOFs, two ZIFs, and four zeolites were studied for removal of SO2, NOx and CO2 from the flue gas mixtures. They used grand canonical Monte Carlo (GCMC) simulations to predict mixture adsorption isotherms and selectivity of the candidate materials for separation of SO2, NOx and CO2 in a mixture containing N2, CO2, O2, SO2, NO2 and NO. In this study, they compared the working capacity, absolute adsorption and adsorption selectivity as three different performance indicators to select the best performing materials. It was concluded that Cu-BTC and MIL-47 were the best adsorbents for separation of SO2 from the flue gas mixture. For the removal of NOx however, Cu-BTC was identified as the best performing material. Finally, for the simultaneous removal of SO2, NOx and CO2, Mg-MOF74 was found to be the best candidate.
The second study from this group was published by Huck et al. [49] in 2014 focusing on screening of 60 different synthesized and hypothetical materials including MOFs, zeolites and porous polymer networks (PPNs) for post-combustion carbon capture. Acknowledging that several evaluation criteria have been already introduced, this publication emphasized the use of parasitic energy as a more realistic metric for materials screening. This is because parasitic energy takes into account the energy penalty associated with the compression process (needed for CO2 storage), as well as several essential thermodynamic properties such as the thermal energy required for heating up the adsorption bed, and the heat required to regenerate it [49]. Using parasitic energy as the evaluation metric, the authors identified Mg-MOF-74, PPN-6-CH2TETA, and PPN-6-CH2DETA as the most promising materials for CCS in coal, natural gas, and direct air capture, respectively.
In a more recent study focused on membrane separation, Qiao et al. [40] screened 137,953 MOFs in an attempt to identify the best performing candidates for separation of CH4, N2 and CO2 using membrane separation technology. In a four-stages strategy, the authors employed a combination of geometric pore characterization metrics (e.g. pore limiting diameter (PLD), percentage of pore size distribution), equilibrium (Henry's constant) and transport properties (diffusivity and permeability) for materials screening showing that the PLD and percentage of pore size distribution are the two key factors governing diffusion and permeation of different gases in MOFs [40].
Almost the entire studies reviewed up to this point had focused on the use of various simple adsorbent metrics for materials screening. They normally include properties such as the pore limiting diameter (PLD), pore size distribution (PSD), Henry's constant of adsorption (KH), adsorption working capacity (WC), selectivity and micropore diffusion. However, in early 2016, Braun et al. [47] published a new study to explore performance of all-silica zeolites for CO2 capture from natural gas where for the first time, inadequacy of some of the above listed adsorbent metrics for materials screening were highlighted. The study suggested that selectivity and working capacity are not necessarily representative of the economic drivers that chemical process designers actually consider [47]. The authors further argued that the use of these metrics can even be deceptive. As an alternative, the authors developed a new metric called separation performance parameter (SPP), which was designed to represent the economic drivers behind CH4/CO2 separation, and applied this metric to explore separation performance and structure-property relationship of tens of thousands of all-silica zeolites recorded in the International Zeolite Association (IZA) database [72] and the Predicted Crystallography Open Database (PCOD) of hypothetical zeolites [73].
The year 2016 also witnessed publication of more sophisticated screening studies including a major contribution from Snurr and co-workers. One of these publications reported on high-throughput screening of MOFs for CO2 capture in the presence of water [41]. The paper focused on the competitive co-adsorption of water as potentially an adverse issue in the deployment of adsorptionbased CO2 capture technologies. Here, the computational screening was conducted to search for MOFs with high CO2/H2O selectivity. The screening workflow consisted of several steps as described below: initially, the framework charges were computed for 5109 MOFs using the extended charge equilibration method (EQeq) [74] which is an approximate, but computationally affordable technique for this purpose. In the next step, the Henry's constants of all MOFs were calculated using the Widom particle insertion. Following this step, the 15 most selective MOFs were identified based on the ratio of Henry's constant for CO2 and H2O. The resulting pre-screened materials were investigated further using more rigorous simulation techniques. For these materials, partial atomic charges were computed using the Repeating Electrostatic Potential Extracted ATomic (REPEAT) method [75], which is obtained from DFT and is more accurate compared to the EQeq technique. Further, GCMC simulations were carried out to calculate the binary and ternary adsorption of CO2/H2O and CO2/H2O/N2 mixtures for the 15 pre-selected MOFs. GCMC simulated-adsorption isotherms were then used to identify MOFs with the highest CO2 selectivity over both water and nitrogen. This study suggests the importance of electrostatic interactions in describing the H2O-MOF interactions. On this basis, the authors suggested that accurate charge calculation methods are required to conduct similar screening studies. They also demonstrated a correlation between small pore sizes and strong binding of CO2 which can limit adsorption of water at high humidity by preventing the formation of water clusters inside these pores [41].
Later in 2017, the same group published a new screening study to explore multivariate metal−organic frameworks (MTV-MOFs) [76]. The authors constructed a new database of ∼10,000 MTV-MOFs with mixed linkers and functional groups. A GCMC-based high-throughput computational screening method was employed to identify the high-performing candidates for CO2 capture. They showed that compared to their parent MOFs, functionalized structures consistently exhibit better CO2/N2 selectivity; and in most cases even CO2 capacity is improved. This work is particularly interesting as it demonstrated that arrangements of mixed linkers containing different functional groups can result in a combinatorial explosion in the number of possible structures which can then be mined to increase structural diversity and surface heterogeneity of materials space. This extended search space may contain candidate materials with higher potential for CO2 capture.
One of the most important developments in 2015 and 2016 was the publication of two studies by Farooq and co-workers [50,77] focused on screening of porous materials for post-combustion carbon capture using a multiscale simulation workflow. In this approach, different types of molecular simulations (i.e. GCMC and MD) were combined with process modelling and optimization of PSA/VSA systems to capture adsorption and transport processes across both microscopic and macroscopic scales. The idea of constructing a multiscale simulation workflow through combining molecular simulations and process optimization for the purpose of materials screening has been originally presented by Hasan et al. [20]. They used this method for cost-effective capture of CO2 using zeolites as adsorbents. Similar multiscale approach was also used by Banu et al. [21] in the context of hydrogen purification using MOFs. Nevertheless it was the above twine study of Farooq and co-workers [50,77] that attracted significant attention among scientific community to the promising capabilities of this new approach for realistic materials screening. In their main screening study, Khurana and Farooq [50] evaluated the performance of 74 real and hypothetical adsorbents in a 4-step VSA process with light product pressurization (LPP). Process optimizations were carried out to minimize overall energy penalty of the process and maximize its productivity while simultaneously meeting the 95% CO2 purity and 90% CO2 recovery criteria for post-combustion carbon capture. As a result of this study, the authors identified several adsorbents with superior performance over 13X zeolite, the current benchmark and the most studied adsorbent for post-combustion carbon capture.
This new development also provided additional evidence that process-based performance indicators such as process productivity, overall energy consumption and product purity do not directly correlate with the intrinsic properties of adsorbent materials [6,50,[77][78][79] that have been widely used by scientists for materials screening over the last decade, as highlighted throughout this historical perspective. Therefore, it is reasonable to conclude that materials screening must be carried out at the process level where separation performance of each material can be evaluated for a given process using the above-mentioned process-level indicators. The multiscale performance-based screening method discussed above addresses several important pitfalls associated with the traditional techniques where materials screening is performed based on intrinsic evaluation metrics: (1) This approach can confirm whether the important CO2 purity−recovery requirement can be met. ISMM, IFMM and HMM classes of evaluation metrics do not take this requirement into account. (2) It can identify the best performance for each adsorbent across a wide range of operating conditions while simultaneously satisfying the purity−recovery constraint. In contrast, adsorbent-based screening methods usually rank materials for a fixed set of operating conditions. (3) The process-level metrics (e.g. energy consumption and productivity) can be directly related to economic drivers of commercialized carbon capture plants (e.g. capital and operation cost). In contrast to this, the intrinsic metrics commonly used in traditional screening approaches, such as working capacity, selectivity, diffusion, etc., cannot be directly correlated with economic considerations of the process.
The above approach for process-based screening of porous materials became particularly important in light of the experimental evidence which support the predictions of the proposed screening platform. In a pilot plant study, Krishnamurthy et al. [80] demonstrated that the 95% CO2 purity and 90% CO2 recovery targets for post-combustion carbon capture can be achieved in experiment using the very same 4-step VSA cycle with light product pressurization which was originally proposed by Khurana and Farooq [50,80]. In a separate study, Perez et al. [81] also verified the ability of multiobjective optimization techniques to guide the design of PSA/VSA processes. In this study, it was shown that purity-recovery Pareto fronts of CO2 as predicted by process modelling of the 4-step VSA-LPP cycle reasonably reproduce the experimental results [81]. These promising observations attracted more attention to the newly proposed process-based materials screening approaches and their combination with molecular simulation techniques, which would allow computational screening of porous materials. Recent examples of these new group of studies are outlined below: In 2018, Farmahini et al. [51] used a similar multiscale platform by combining GCMC simulation with process modelling and optimization of the 4-step VSA-LPP cycle to explore the challenges associated with the interface between molecular and process levels of description. In this study, they identified several sources of inconsistency for accurate implementation of the multiscale screening workflow which can largely affect correct prediction of materials performance at the process level. This includes the numerical procedures adopted to feed the equilibrium adsorption data into the process simulation, and the role of structural characteristics of adsorbent pellets including pellet porosity and pellet size.
In 2019, Balashankar and Rajendran [82] employed a two-stage approach to screen 119,661 hypothetical zeolites, 1031 zeolite immidazolate frameworks, and 156 zeolites identified by the International Zeolite Association. In their study, the first stage was dedicated to the rapid screening of all materials under investigation using a computationally inexpensive batch adsorber analogue model to filter adsorbents which can meet 95% CO2 purity and 90% CO2 recovery targets. This stage was then followed by detailed process modelling of 15 top-performing candidates from the previous stage in addition to 24 synthesizable zeolites using the widely used 4-step VSA-LPP cycle to estimate the process level performance indicators more accurately. Out of the 39 adsorbents screened in the second stage, 16 material candidates outperformed Zeolite 13X both in terms of productivity and energy consumption [82].
In 2020, a series of new studies were published focusing on development of more sophisticated simulation techniques for multiscale screening of porous materials. Following their previous study from 2018 [51], Farmahini et al. [83] published a new study in 2020 proposing a new approach for estimation of pellet morphology (e.g. pellet size and porosity) that cannot be calculated from molecular simulations as part of the multiscale screening workflow. This study is particularly important for consistent screening of porous materials, considering variation of pellet morphology can greatly alter separation performance at the process level. The authors demonstrated that a series of competing mechanisms associated with diffusion into adsorbent pellets, convection mass transfer through adsorption column, and pressure drop across the bed can be tuned through optimization of pellet size and pellet porosity to maximize separation performance of different classes of porous materials including zeolites and MOFs [83].
Later in the same year, Burns et al. [84] screened 1,632 experimentally characterized MOFs using a multiscale platform which combines molecular simulations with process optimization and machine learning models. In their screening study, they also employed the famous 4-step VSA-LPP cycle and found that a total of 482 materials can achieve the 95% CO2 purity and 90% CO2 recovery targets out of which 365 materials have parasitic energies below that of commercial solvent-based CO2 capture technologies [84].
Another screening study in 2020 was published by Yancy-Caballero et al. [85] who compared process level performance of 15 promising MOFs with Zeolite 13X as a benchmark using three different process configurations including a modified Skarstrom cycle, a five-step cycle, and a fractionated vacuum swing adsorption cycle. The results from this study suggests that UTSA-16 and Cu-TDPAT perform equally well or even better than Zeolite 13X in all the three process configurations mentioned above. The authors also compared process-level ranking of these MOFs with other rankings obtained based on simplified HMM and GEM metrics. They showed that the rankings suggested by these metrics may differ significantly from the one predicted by detailed process optimization [85].
In this context, the most recent study in 2020 was published by Pai et al. [86] who developed a generalized and data-driven surrogate model which can reproduce operation of PSA/VSA processes at cyclic steady state with high accuracy. The multiscale screening framework developed here simultaneously optimizes adsorption isotherm properties and process operating conditions in order to estimate performance indicators of the process. The framework makes use of a dense feed forward neural network trained with a Bayesian regularization technique and is able to significantly reduce the simulation and optimization time required for multiscale screening of porous materials for postcombustion carbon capture [86].

Multiscale Screening Workflow
In the previous sections, we briefly discussed why materials screening is important in the context of PSA/VSA technologies for the post-combustion carbon capture. We also provided a historical perspective on the evolution of materials performance metrics and screening methods, which have been used so far. The main objective of these sections was to illustrate to the reader the importance and the gradual evolution of the research community towards adopting more complex multiscale screening workflows as the emerging way to realistically evaluate separation performance of porous materials. The objective of the current section however, is to introduce the essential elements and methods involved in the multiscale screening workflows. The structure of this section logically follows the multiscale workflow diagram, shown in Figure 6. The starting point of this workflow is a database of porous materials. In Section 6.1, we review the currently existing materials databases and the computational tools required to characterise structural properties of the porous materials from these databases.
Molecular simulations are used to obtain equilibrium and transport properties at a molecular level. These methods are introduced in Section 6.2. Finally, following the workflow we pass the information from molecular simulations to the process level modelling and optimization. Models, methods, and data required for this stage are reviewed in Section 6.3.

Material Databases and Characterization Tools
This section corresponds to the first step in the multiscale material screening workflow. The aim here is to provide concise and practical reference to the reader on what databases are currently available, what materials and data they contain and what tools are available to build geometric descriptors for materials in these databases.

Databases of Porous Materials
MOFs are the primary and most prominent example of the emerging families of materials and it is useful to briefly review what these materials are. Although the origins of MOFs can be traced as far back as the late fifties, they were given their current name, Metal-Organic Frameworks, in the seminal paper by Yaghi and Li in 1995 [87]. To prepare a MOF one uses two types of building blocks: metal complexes and organic molecules capable of forming strong coordination bonds with these complexes. In the synthesis process, the building blocks form a crystalline framework where metal complexes comprise the vertices of the framework, connected by the organic linkers. Several papers that followed in the late nineties discovered few more examples of these frameworks, however, most importantly they demonstrated that these structures possessed permanent stable porosity, high surface area and new materials could be designed simply by variation of the building blocks, leading to the concept of isoreticular material design [88][89][90]. Since then, tens of thousands of new MOFs have been discovered: the most current assessment of the Cambridge Structural Database suggests ca. 100,000 reported structures that can be qualified as MOFs [91], while the modular nature of MOFs implies that in principle infinite variation of structures is possible (if we assume that the diversity of MOFs can approach the diversity of the organic chemical space).
ZIFs, discovered a few years later [92,93], is a subclass of MOF materials that have zeolite framework topologies in which silicon atoms are replaced by transition metals and the bridging oxygens are substituted by imidazolate building units [94]. Currently, there are about 300 ZIFs reported in the CSD and potential application of these materials in the context of chemical separations has been recently reviewed by Pimentel et al. [95]. In contrast to materials based on coordinative assembly and coordination bonds, Covalent Organic Frameworks (COFs) do not feature metal complexes and are based on covalent bonds [96]. Since their discovery in 2005, a substantial number of 2D and 3D COFs have been reported with diverse structural and chemical properties [96].
Crystalline materials, such as MOFs, ZIFs, COFs can be contrasted with several traditional and emerging classes of amorphous porous polymers, such as activated carbons [97], carbide-derived carbons [98] and Polymers with Intrinsic Microporosity (PIMs) [99]. Porous Aromatic Frameworks (PAFs) is another class of porous materials with rigid aromatic open-framework structure constructed by covalent bonds. Although, PAFs are not crystalline they are ordered with regular and high porosity [100].
This wealth of new materials should not overshadow more traditional classes of porous materials such as zeolites, which, due to their stability, attractive cost, commercial availability and maturity in industrial applications, will likely remain the primary adsorptive materials for years to come. There are currently ~200 zeolite topologies recognized by the International Zeolite Association. Using computational methods, millions of new topologies have also been hypothesized [73,101]. An ongoing research is to understand the magnitude and diversity of the materials landscape for adsorption science [102,103] and to evaluate what portion of this structural space is realizable in experiments [104]. Combined, these classes of materials provide an enormous chemical and structural diversity, collectively described as the materials genome [105].  [111,112] and Gómez-Gualdrón et al. [108] generated two new databases of hypothetical MOFs using topology-based algorithms that are topologically more diverse. These two databases are called BWdatabase and ToBaCCo databases throughout this paper.

(a) Wilmer's Database:
This database contains 137,953 structures and is generated by recombining a library of 102 building blocks including secondary building units (SBU) and organic linkers from crystallographic data of already synthesized MOFs using a "Tinker-toy" algorithm [61]. The resulting hypothetical database is however composed of a few underlying framework topologies [113] which is due to the use of building blocks that are topologically largely similar [113]. By testing a very limited set of MOFs including HKUST-1, IRMOF-1, PCN-14 and MIL-47 , the authors suggested that their method can closely reconstruct atomic structures of the experimentally synthesized materials [61]. Nevertheless, generalization of this finding is subject to more comprehensive validations, considering no energy minimization has been performed for any of the constructed structures in this database. Also, the materials predicted in this database do not include electrostatic charges, hence their applications will be limited to very few adsorption cases where electrostatic interactions are not important (e.g. CH4 adsorption). In fact, Wimer et al. used this database to search for MOF materials suitable for methane storage. They identified more than 300 MOFs with a predicted methane-storage capacity larger than that of any previously known material [61].

(b) BW-Database:
This new database of hypothetical MOFs was constructed using the topology-based algorithm of Boyd and Woo [111] and contains 324,426 structures which are generated by assembling a set of secondary building units containing 8 inorganic and 94 organic SBUs resulting in 12 different topologies [112]. The set was further diversified by chemical modification of MOFs, in which available hydrogens were replaced by functional groups. All MOFs in this database are structurally optimized using classical force fields. Framework charges for all structures included in this database were also computed using the charge equilibration method (QEq) and the MEPO parameters [112].
(c) ToBaCCo Database: This database constructed using the Topologically Based Crystal Constructor (ToBaCCo) algorithm and contains 13,512 MOF structures with 41 different edge-transitive topologies [108,114]. The database makes use of a top-down construction algorithm which uses topological blueprints and molecular building blocks as input to assemble MOF structures. The algorithm does not check for atom overlaps as part of the construction process therefore the resulting structures must be structurally optimized before being used in molecular simulations [114]. It also does include partial atomic charges.

Cambridge Structural Database (CSD):
The Cambridge Structural Database (CSD) contains more than a million organic and metal-organic small-molecule crystal structures which are obtained from X-ray or neutron diffraction analyses [106]. The MOF structures deposited in this database are experimentally realized, nevertheless the use of CSD entries for high-throughput screening of porous materials is not straightforward. Checks must be performed to make sure that the candidate structures obtained from CSD are adequately porous and are free from residual substances that are leftover of the synthesis processes. As such, the first step in performing high-throughput screening of experimental MOFs is to construct curated subsets of CSD that can fulfil the above criteria.

Goldsmith Database of Experimental MOFs:
In 2013, Goldsmith et al. [115] constructed a MOF database containing 22,700 computation-ready structures which was derived from the CSD after the removal of unbounded guest molecules (e.g. residual solvents). By excluding disorder compounds and those with missing atoms, the total number of MOF structures were reduced to 4,000 [115] which do not include those with interpenetrated frameworks and charge-balancing ions [107]. The materials included in the database were subsequently characterized by calculating porosity, surface area and total theoretical H2 uptake [115]. Goldsmith et al. used their MOF database to estimate the maximum theoretical uptake of hydrogen based on the so-called "Chahine rule" [116] known for hydrogen adsorption in microporous carbons but also shown to be valid across a wide range of other porous materials including MOFs [115].

CoRE-MOF Database: Construction of the Computation-Ready, Experimental Metal-Organic
Frameworks (CoRE-MOF) was a major attempt in development of a MOF database that can be directly used in molecular simulations. The first version of CoRE-MOF [107] contains 5,109 3D MOF structures with pore-limiting diameter greater than 2.4 Å which are derived from CSD. The MOF structures were screened to make sure that all MOFs included in the database are crystalline (no disorder) and solventfree. The database also reports helium void fractions of all MOFs in addition to their surface area, accessible volume, largest cavity diameter (LCD) and pore-limiting diameter (PLD). In the original version of the database, the structures were not optimized (except for very few MOFs that were manually edited) [107]. Following the initial release of CoRE-MOF, two modified subsets of this database were released in 2016 and 2017. The first subset contains 2,932 experimental MOFs whose partial atomic point charges were calculated using planewave DFT and the DDEC charge partitioning methods [117]. The second subset focuses on the geometry optimization of 879 experimentally synthesized MOFs using a periodic density functional theory (DFT) method [118]. The latter publication demonstrated that although the majority of MOF structures undergo less than 10% change in their structural parameters (e.g. pore size, lattice parameters, unit cell volume, helium void fraction) upon DFT optimization, many other MOF structures change significantly after geometry optimization especially those materials whose crystalline structure were cleaned from solvent residues. More importantly, it was shown that the DFT optimization had a large impact on simulated gas adsorption in some cases, even for materials whose crystalline structure did not change significantly [118]. This study has important implication for high-throughput materials screening approaches that rely on databases of experimentally synthesized materials such as CSD [106]

CSD-MOF Subset:
In 2017, Moghadam et al. [110] constructed a new subset of CSD for solvent-free MOFs in which 69,666 1D, 2D and 3D MOFs were listed out of which 54,808 structures are nondisordered. These materials were characterized using the Zeo++ code [63] based on the Voronoi decomposition technique to calculate the accessible surface area, accessible pore volume, LCD, and PLD. It was found that 46,420 structures have gravimetric surface area equal to zero which essentially means that N2 size molecular probes cannot access their pore spaces for geometric surface area calculations [110]. It is shown that the remaining 8,388 MOFs have PLD values larger than 3.7 Å which is approximately 3,600 structures more than what was previously published by Chung et al. [107] in the initial version of the CoRE-MOF database. Currently, the MOF subset of CSD database contains approximately 100,000 MOFs [91]. The main advantage of the CSD-MOF subset is that it is integrated into the Cambridge Crystallographic Data Centre's (CCDC) structure search program. This not only allows for tailored structural queries (e.g. generation of MOF subsets based on secondary building units or selection of non-disordered materials), it can also be used to automatically update the database with subsequent addition of new MOFs to CSD [110].

Hypothetical Zeolites Database (hZeo-DB):
hZeo is a database of computationally predicted zeolitelike structures which are generated by systematically exploring 230 space groups, unit cell dimensions between 3 Å to 30 Å, and T-atom densities from 10 to 20 per 1000 Å 3 [73,120]. A computational procedure based on Monte Carlo search was employed to produce 3.3 million zeolite-like structures out of which 2.6 million topologically distinct structures were identified after energy minimization [120]. Roughly 10% of this number are the structures deemed to be thermodynamically accessible as aluminosilicates based on energy stability of the structure [73].

Database of Zeolite Structures (IZA-DB):
IZA-DB provides information about the structure of all the zeolite framework types that have been approved by the Structure Commission of the International Zeolite Association (IZA-SC). The database currently contains 241 ordered and 11 partially disordered topologies [72].

Database of Hypothetical Porous Polymer Networks (hPPN-DB):
The hypothetical PPN databse constructed by Martin et al. [109] contains almost 18,000 hypothetical structures of porous polymer networks which are predicted in silico using commercially available chemical fragments and two experimentally known synthetic routes; hence aiming to provide a database of synthetically realistic PPNs [109]. All structures from this database have their structure optimized using semiempirical electronic structure methods [109]. The structures are also characterized for their topological properties (e.g. accessible surface area, accessible volume, diameter of the largest included sphere) and methane adsorption characteristics [109].

Hypothetical COF Database (hCOF-DB):
This database is a collection of 69,840 novel hypothetical covalent organic frameworks (COFs) which are assembled from 666 distinct organic linkers and four established synthetic routes [121]. It contains 18,813 interpenetrated 3D structures, 42,386 noninterpenetrated 3D structures and 8,641 2D-layered structures. All materials are structurally relaxed using classical force fields. The database does not include partial atomic charges for the deposited COFs.

CoRE-COF Database:
In 2017, Tong et al. [122] compiled a computation-ready database of experimental covalent organic frameworks (COFs) containing 187 structures. The original version of the database contains 19 3D-COFs and 168 2D-COFs. The structures collected in this database are reported to be disorder-free and solvent-free which make them ready for computational studies. Although most of the structures available in CoRE-COF database are cleaned versions of experimentally reported CIF files, some of the COFs collected here are constructed computationally based on the information reported in the literature where synthesis of the corresponding COFs had been reported without any CIF file. CoRE-COF materials are structurally optimized and refined using a two-step procedure with classical force fields and the dispersion-corrected DFT method of Grimme (DFT-D2) [122]. The database also reports on structural features of each COF including their largest cavity diameter, pore-limiting diameter, accessible surface area and free volume. Since its first release, the CoRE-COF database has been updated regularly so that its most recent version (CoRE-COF Ver. 4.0) 1 contains 449 structures with framework charges included as obtained from the charge equilibration (Qeq) method.
CURATED COF Database: Clean, Uniform, Refined with Automatic Tracking from Experimental Database (CURATED) of covalent organic frameworks (COFs) is another database of experimentally realized COFs. The initial version of the database includes 324 structures, however the database has been updated recently so that its most recent version (Feb 2020) contains 482 structures. All structures collected in CURATED COFs are cleaned from solvent molecules and have no partial occupation or structural disorder. They are structurally optimized using DFT with DDEC framework partial charges included [123].

Hypothetical ZIFs Database (hZIF-DB):
In 2012, Lin et al. [48] published a paper on computational screening of large number of zeolites and zeolitic imidazolate frameworks (ZIFs) for carbon capture.
In this study, ZIF structures were generated computationally by using zeolite topologies of the International Zeolite Association (IZA) database. In doing so, the distance between zinc atoms and the centre of imidazolate rings was set to be 1.95 times larger than the silicon-oxygen distance in zeolites. ZIF frameworks were then generated by scaling the corresponding zeolite structures by the same factor and replacing every oxygen atom with an imidazolate group and substituting every silicon atom with a zinc atom. The resulting ZIF geometries were validated by comparing against geometries of two experimentally known ZIF structures (i.e. ZIF-3 and ZIF-10) [48]. This structure is not, however, available online or in a depository to further comment on its characteristics.

Database of Porous Rigid Amorphous Materials (PRAM-DB):
Hitherto, there have been many efforts for systematic development of various materials databases for crystalline porous solids. Nevertheless, attempts for construction of such databases for amorphous porous materials have been scarce. In an important development, Thyagarajan and Sholl [126] have recently collected 205 atomistic models of amorphous nanoporous materials which had been previously published by various groups. This new database of porous rigid amorphous materials (PRAM-DB) contains several classes of materials with disordered porous structures including amorphous zeolite imidazolate frameworks (a-ZIFs) [127], activated carbons [128], carbide-derived carbons [129][130][131][132][133], polymers of intrinsic microporosity (PIMs) [134][135][136][137], hyper-cross-linked polymers (HCPs) [138][139][140], kerogens [141] and cement [142] which all have important applications in adsorption separation technologies. The database contains partial atomic charges for most of the materials. It also reports on a wide range of physical properties for each material. This includes pore limiting diameter (PLD), the largest cavity diameter (LCD), the accessible surface area and pore volume, pore size distribution (PSD), ray-tracing histograms, PXRD patterns and radial pair distribution functions (RDF) [126]. The new study also reports single-component and binary adsorption isotherms of several gases for these materials [126]. As can be seen from the reviewed studies, classification of materials within the databases and early efforts in computational screenings are based on the geometric descriptors of porous materials, such as the accessible surface area, pore limiting diameter and pore volume. As this is a practice-oriented review, we believe it is useful to mention the material characterization software available to obtain these geometric properties for crystalline and amorphous porous structures. To begin with, we refer the reader to several articles describing what properties of porous materials can be calculated computationally and how they are related to the properties that can be measured and to the physical process of adsorption in porous materials [63,64,143,144]. In principle, calculation of selected properties, such as the solvent-accessible surface area (in application to porous materials often called simply, accessible surface area), is available within many commercial and free software packages. Three packages available for a more comprehensive assessment of the materials are Poreblazer [64,143,145], Zeo++ [63], and PorosityPlus [146]. From this list, Poreblazer developed by Sarkisov and Harrison [64,145] and PorosityPlus developed by Opletal [146] are written in Fortran and are available as open-source packages. Zeo++, developed by Haranczyk and co-workers [63] is a C++ package based on the Voronoi tessellation methods [63]. With Voronoi network being a dual graph of Delaunay network, the approach employed by Zeo++ is closely related to that of Foster et al. [147]. The program is downloadable from the website of the developers, with the source code available upon request only.
All three codes mentioned above are able to calculate accessible surface area (equivalent to the area of the surface formed by the nitrogen probe rolling on the surface of the atoms of the structure), pore volume (using several alternative definitions of this property) and pore size distribution. Poreblazer and Zeo++ can also calculate pore limiting diameter (PLD) of the porous frameworks, while PorosityPlus is able to compute radial distribution function (RDF) of the adsorbed phase in the system. One important feature of Zeo++ software is its ability in reading framework structures in CIF format, while the other two programs can only use XYZ format as their input for the porous framework. A detailed comparison of Poreblazer, Zeo++ and RASPA [148] has been recently provided by Sarkisov et al. [145] for structural characterization of CSD-MOF Subset database [110]. Here, we note that RASPA is a molecular simulation software which is mainly known for its capabilities for Monte Carlo simulations. This program is presented in the following section where we discuss grand canonical Monte Carlo (GCMC) technique for simulation of equilibrium adsorption isotherms.

Molecular Simulation
The purpose of this section is to briefly introduce the two main and most widely used molecular simulation techniques, grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulations which are respectively used for simulation of adsorption and transport properties on a microscopic level. In particular, we would like our intended reader to appreciate what parameters are required for these simulations and their sources, open-source software available to perform these simulations, as well as current challenges and gaps in the field.
In the last 50 years molecular simulations have become an important tool of enquiry in chemical engineering, material science, chemistry and other fields. This is driven by the development of advanced simulation codes and the availability of computational power. Furthermore, concurrent advances in computational quantum chemistry made it possible to consider complex, multiatomic systems at a quantum level, and in this way develop more accurate force fields for representation of the system in classical mechanical simulations. In the context of adsorption problems, comprehensive reviews on molecular simulations for Metal-Organic Frameworks have been provided by Yang and coworkers [149] and for zeolites by Smit and Maesen [150].
In Section 6.2.1, we introduce fundamentals of GCMC method followed by Section 6.2.2 which is dedicated to introducing the main publically available simulation software for performing this type of simulation. Next, in Section 6.2.3, fundamentals of molecular dynamics will be discussed, which will be also followed by a section related to the open-source programs that can be used to run MD simulations (Section 6.2.4). Finally in Section 6.2.5, we will briefly introduce atomic force fields which are central to accurate simulation of molecular systems and will review the current gaps in the field of force field development and comment on their implications for multiscale materials screening studies.

Fundamentals of Grand Canonical Monte Carlo Simulations
In this section, we briefly review the grand canonical Monte Carlo simulation method, which is widely used for calculation of equilibrium adsorption data. For a more comprehensive review of Monte Carlo methods, we would refer the reader to reference books [151][152][153] and several excellent articles by Dubbeldam and co-workers on the Monte Carlo methods and the organization of computer codes associated with them [148,154].
The problem of interest here is the adsorption of small molecules (carbon dioxide, nitrogen, methane, hydrogen) in crystalline porous materials. The volume, V , and temperature, T, of the system are fixed, and the specified value of the chemical potential, μ, establishes thermodynamic equilibrium between the system and the bulk reservoir, serving as a source and sink of adsorbate molecules. From the statistical-mechanical point of view, the system corresponds to the grand-canonical ensemble (μ V T), for which the Metropolis Monte Carlo method serves as a conventional simulation technique of choice. This approach is suitable for rigid porous materials, which do not exhibit significant changes in the specific volume upon adsorption and desorption of guest molecules. For flexible materials, which represent a significant and interesting class of structures within the MOF families, more advanced simulation methods exist (such as Osmotic Ensemble Monte Carlo and Gibbs Ensemble Monte Carlo [154]).
Within the GCMC method, configurations of the system are generated via a set of standard trial moves; translation, rotation (in case of rigid molecular species), insertion and deletion, with the following acceptance probability applied to ensure the Boltzmann distribution of the generated states: a) Translation: where U represents the potential energy, , and V are the number of molecules and volume respectively, is the reciprocal thermodynamic temperature, 1/kBT, with kB being the Boltzmann constant; is an Euler angle of the rigid body rotation, f is the fugacity of the adsorbing species, which is related to the chemical potential as: where is the rotational partition function for a single rigid molecule, equal to 1 for a single particle molecule, and Λ is the thermal de Broglie wavelength: where ℎ is the Planck's constant and m is the molecule mass. The schematic diagram of the GCMC workflow is shown in Figure 7. According to this scheme, a Monte Carlo simulation of adsorption requires the following inputs: -Force field parameters: these parameters define what atoms and molecules are present in the system and describe how they interact with one another. This includes parameters associated with non-bonded dispersion interactions, partial charges on the atoms of the structure and molecules, geometry of the adsorbing molecules (distances and relative positions of the atoms within the molecule).
-Initial configurations of the species present in the system: this includes positions of the atoms of the porous crystal structure and positions of any already adsorbed molecules.
-Simulation parameters, which specify various aspects of the actual simulation run, including the number of Monte Carlo moves and their weights; number of steps allocated for the equilibration of the system; parameters associated with the statistical analysis of the simulation (i.e. number and size of blocks in the block-average analysis), temperature conditions and fugacities of the adsorbing components. This input data category may also prescribe particular specialized methods to calculate electrostatic interactions between partial charges on the atoms.
These data are submitted to the Monte Carlo simulation engine, where the microstates of the system are stochastically changed via one of the available Monte Carlo moves. The probability to accept the move is calculated according to equations (1) -(4). As the simulation progresses, the positions of the molecules change and the number of the molecules fluctuates, producing an ensemble of microstates over which the average properties of the system can be calculated. This ensemble of the microstates is called a trajectory and it is a common outcome of both Monte Carlo and molecular dynamic simulations (in a sense that it reflects the position of the system in the phase space), with the difference that the Monte Carlo trajectory is not continuous and does not contain information about the velocities of the molecules.
The ensemble of microstates within the trajectory can be used to produce the relevant output properties of the system. In the context of the adsorption studies, the most important property is the average number of molecules present in the system, which when plotted as a function of pressure values produces an adsorption isotherm. An important distinction has to be made between the absolute and excess amount adsorbed. The absolute amount adsorbed is the actual number of molecules present in the micropores at a particular fugacity. The excess amount is the difference between the absolute amount adsorbed and the number of molecules that would be present in the micropore volume according to the bulk gas density at the pressure and temperature of adsorption. The distinction between different definitions of adsorption and their connection to the experimental measurements has been discussed by Brandani et al. [155]. Monte Carlo simulations report absolute mount adsorbed, whereas experimental measurements are more often presented as excess amount. The process simulations discussed in the next section take as an input analytical models for the absolute amount adsorbed.
Another important property that can be obtained from GCMC simulation is enthalpy of adsorption. As will be discussed in the process modelling section of this review, in real processes and in process models based on adiabatic considerations, heat effects may play a role in the performance of the cycle. In molecular simulations, this property can be calculated either using the expression based on the result from the statistical-mechanical fluctuation theorem [156] or, in a direct analogy to the experimental methods, using the Clausius-Clapeyron equation. In the first case, a single isotherm is sufficient to calculate the heat of adsorption at each adsorption pressure. However, at high loading the reliability of this method deteriorates: it relies on the fluctuation of the number of adsorbed molecules in the system, and since at high loading the acceptance ratio for the insertion and deletion Monte Carlo moves is low, convergence of the method becomes problematic. This is not an issue for the approach based on the Clausius-Clapeyron equation [156], however, this method requires adsorption isotherms at several temperatures. Finally, simplified expressions are available if one is interested in the heat of adsorption in the Henry's law (zero loading) regime.
In addition to the properties directly required by the process simulation data (adsorption equilibria, heats of adsorption), molecular simulations also generate a wealth of information by visualizing the adsorption process on a molecular level (e.g. visualizations and density maps). These properties help to elucidate, for example, presence of specific binding sites and distribution of the molecules in the structure, which in turn can be used to construct new analytical models for adsorption.
So far, this brief introduction to the grand canonical Monte Carlo methods for adsorption problems implicitly assumed rigid crystal structures and rigid adsorbate molecules (with small gas molecules, such as nitrogen, carbon dioxide, methane being adequately described by this approximation). Extension of GCMC simulations to larger flexible molecules (i.e. alkanes) requires more advanced techniques, such as the Configurational-Bias GCMC [150]. Adsorption behaviour in flexible MOFs has attracted significant attention over the years. To capture these phenomena, simulation in the osmotic ensemble is required as well as advanced force fields to correctly represent the internal degrees of freedom within the framework [157].

Monte Carlo Simulation Codes
To make the review a practical reference, here we briefly introduce the open-source Monte Carlo codes for simulation of equilibrium adsorption isotherms in porous materials. These codes are listed in Table 4. We note here that a special issue of Molecular Simulation journal invited the community to reflect on the codes and algorithms available for the Monte Carlo simulations, their accessibility and applicability, efficiency and challenges [158]. In a recent study, we tasked ourselves with exploring the consistency of some of the most commonly used MC codes as listed in Table 4 and examined their relative efficiency [159]. For this, we concentrated on a specific case study of carbon dioxide adsorption in IRMOF-1 material at conditions for which previous simulation results and experimental data were available [160]. It was a significant reassurance for us to observe that the codes are indeed consistent with each other. To assess their relative efficiency, we employed analysis based on the statistical inefficiency of sampling to compare trajectories from different codes on a consistent basis of the rate with which they were generating a statistically novel configuration. Our analyses revealed some differences in the overall performance of various MC codes, nevertheless this variation was found to be relatively negligible [159]. RASPA, MuSiC and DL_MONTE were overall the top performing programs in the analysis. Within the same article, we also generated consistent setups and scripts for all the codes for the above test case, which can be used by the molecular simulation community as a template for consistency tests and validation of future MC codes. These materials are available free of charge from our online Bitbucket repository 2 [159]. Consistency and efficiency of MC codes are particularly important in the context of materials screening and multiscale simulation workflows.  [161] https://cassandra.nd.edu/ DL Monte Purton et al. [162] https://www.ccp5.ac.uk/DL_MONTE MuSiC Gupta et al. [163] http://zeolites.cqe.northwestern.edu/Music/ RASPA Dubbeldam et al. [148] https://www.iraspa.org/RASPA/index.html Towhee Martin [164] http://towhee.sourceforge.net/ Here, we briefly introduce the codes listed in Table 4: Cassandra: is a MC program developed in the group of Maginn at the University of Notre Dame. It is an effective package for simulation of the thermodynamic properties of fluids and solids [161]. Cassandra supports canonical (NVT), isothermal-isobaric (NPT), grand canonical (μVT), osmotic (μpT), Gibbs (NVT and NPT versions) and reactive (RxMC) ensembles. The code can be compiled to run in parallel using OpenMP [161].
DL_Monte: is another Monte Carlo simulation software that can be run in parallel [162]. It is originally developed by Purton and co-workers at Daresbury Laboratory in the UK with special emphasis at materials science. It is now being developed as a multi-purpose simulation package in collaboration with Wilding (University of Bristol) and Parker (University of Bath) research groups. The code can simulate systems in canonical (NVT), isobaric-isothermal (NPT), grand canonical (µVT), semi-Grand canonical, and Gibbs ensembles [162]. DL_MONTE is a twin sister code of DL_POLY package, a molecular dynamics simulation software that will be introduced later in this review. With regard to parallelization of MC codes such as DL_Monte and Cassandra, Gowers et al. [159] have demonstrated that the measured performances of existing implementations show poor efficiency due to various reasons. At least in the context of adsorption simulations and computational screening of porous materials, parallel execution of multiple MC runs offers higher efficiency and larger overall speed up as compared to parallelization of MC codes [159].

MuSiC:
The Multipurpose Simulation Code (MuSiC) is an object-oriented software developed in Snurr's research group from Northwestern University in the US [163]. The code supports grand canonical (μVT), canonical (NVT), and isobaric-isothermal (NPT) ensembles. It can also be used to perform hybrid MC and molecular dynamics (MD) simulations [165].
RASPA: is a molecular simulation code designed for simulation of adsorption and diffusion processes in flexible nanoporous materials [148]. The code was developed in Snurr's group at the Northwestern University in collaboration with several other scientists in the field of molecular simulations [148]. RASPA supports a variety of ensembles including micro canonical (NVE), canonical (NVT), isobaricisothermal (NPT), isoenthalpic-isobaric (NPH), Gibbs (NVT and NPT versions), and isobaric-isothermal ensemble with a fully flexible simulation cell (NPTPR) [148]. It can be used to perform both Monte Carlo and molecular dynamics simulations, however it is best known for its capability as a MC software. The code also supports configurational bias Monte Carlo (CBMC), and continuous fractional component Monte Carlo (CFMC) for rigid and flexible systems [148,154].

MCCCS Towhee:
The Monte Carlo for Complex Chemical Systems (MCCCS) program was originally developed in Siepmann's research group at the University of Minnesota. It is currently being developed and maintained by Martin [164,166]. The code was initially designed for the prediction of fluid-phase equilibria, however it has been extended later to simulate different systems including porous materials. Towhee supports a variety of ensembles including NVT, NPT, μVT and Gibbs [164].

Fundamentals of Molecular Dynamic Simulations
In this section, we turn our focus to molecular dynamics, which is widely employed for calculation of time-dependent phenomena across different fields from gas separation to materials science, geological sequestration of gases, biomolecular science, and drug discovery [149,[167][168][169][170][171][172]. The brief description provided here solely concerns molecular diffusion of simple gases in crystalline porous materials, which is relevant to the topic of this review. This section is meant to serve as an introductory material for non-expert readers. For more in-depth discussion of this technique, the reader is referred to numerous resources available in the literature [151,[172][173][174][175][176].
Accurate and quantitative description of molecular transport in porous materials is crucial for the design and development of separation processes. This is normally attempted using the transport theories that describe the mechanisms of mass transport under confined pore spaces [177,178]. Molecular dynamics however provides an alternative route for calculation of gas transport inside microporous materials without dealing with the uncertainties arising from structural heterogeneity of real samples in experiment [177]. MD can be used where classical theories fail [175,177,179] or to check the validity of analytical models when experimental data are hard to obtain [180][181][182][183]. It is a computational technique for calculation of the equilibrium and transport properties of classical manybody systems [173]. In MD, we choose a system consisting of N particles for which we continuously solve Newton's laws of motion at constant time steps to explore time evolution of the potential energy surface. Therefore, in contrast to Monte Carlo method where molecular movements of the system are modelled stochastically, in MD molecular motions and hence their time-dependent positions in space are precisely calculated based on Newtonian mechanics [175]. Here, the total force exerted on each particle is given by [175,184] ( ) = ( ) = = 2 2 = − ( ) = 1, 2, 3, … , where , , , and denote the force, velocity, mass, acceleration and position of the i th particle respectively. and stand for total potential energy and time. The above equation is normally solved from a Taylor series expansion about initial position and velocity of particles in the system [175,185].
There are several algorithms in the literature for time integration of equation (7) such as the Leapfrog [186] and Verlet [184]. In the latter one, which is not only one of the simplest methods but also one of the most widespread algorithms [173], the position of the particles at each time step is calculated by: The above estimate of the new position of particle i contains an error that is of order ∆ 4 , where ∆ is the time step in the MD simulations [173].
In the context of gas adsorption where diffusion of particles in porous materials is monitored, MD simulations are normally carried out in the canonical (NVT) ensemble where volume, V , temperature, T, and the number of particles in the system, N, are conserved. This approach is suitable for molecular diffusion in materials whose porous frameworks are rigid. For materials with large framework flexibility, simulations can also be performed in NPT ensemble where pressure, P, is constant instead of the system volume, V [187,188]. This would allow volume of the system to change under constant pressure, which is often the case in diffusion experiments. Figure 8 depicts the schematic diagram of the MD workflow and the properties that can be calculated from typical MD simulations. In MD, we need to define a set of starting (i.e. initial) configurations for the system which are often obtained from GCMC simulation. Similar to the MC method, interatomic interactions of all particles must be defined using an appropriate set of force fields along with other simulation parameters that are normally supplied to an MD program as input data (e.g. time step, temperature, pressure, etc). MD generates time trajectory of the system containing atomic positions of all particles and their associated potential energies. From these data, a number of transport [150,172,189], and thermodynamic properties [167,[190][191][192] can be calculated. In this section however, we only focus on transport properties (i.e. molecular diffusivities) that can be computed using MD.
Diffusion of adsorbate molecules through the porous structure of adsorbent materials can be conveniently calculated using MD simulations, however measurement of experimental data for molecular transport is a challenging task particularly for mixtures [176,193]. Here we introduce several key concepts associated with the diffusion of molecules in porous media and molecular dynamics as a technique to study this processes, first considering a single component case and then exploring the formalism for the multicomponent systems.

Diffusion in single-component systems:
For single-component systems, self-diffusivity, collective diffusivity, and transport diffusivity are three types of diffusion phenomena that are commonly studied by molecular simulations [150,172,176,189,[194][195][196][197][198]. Self-diffusivity (Ds) describes the motion of individual labelled molecules through a fluid in the absence of the chemical potential or concentration gradients. In experiments, this property is measured using tracer diffusivity techniques, such as pulsed field gradient (PFG) NMR. In equilibrium molecular dynamics (EMD), self-diffusivity can be calculated from the mean-squared displacement of labelled particles using the Einstein relationship given by: Where d is dimensionality of the system. Ds can also be computed from the time integral of the velocity autocorrelation function (VACF) defined by: here, vi (t) is the centre of mass velocity vector of molecule i. The brackets in (9) and (10) indicate an ensemble average taken over the simulation run time. As diffusion in porous materials is an activated process, temperature dependence of is typically captured in the well-known Arrhenius relation where 0 is the pre-exponential constant and is the activation energy.
In contrast to self-diffusivity, the transport (Dt) and collective, or corrected, (Dc) diffusivities are associated with the macroscopic flux of molecules arising from the spatial concentration gradient in the fluid [149,175]. The transport diffusivity also referred to as the Fickian or chemical diffusivity is related to net flux in the system, which is described by Fick's first law: here, J and ∇q are the flux and concentration gradient in the adsorbed phase, respectively. Equation (11) can also be described in terms of the chemical potential gradient ∇μ [175]: where, L is the Onsager transport coefficient, and is the collective diffusivity which is also known as the corrected or Maxwell-Stefan diffusivity [175].
The transport diffusivity ( ) is related to the collective diffusivity ( ) through a term associated with curvature in the adsorption isotherm [189]. This parameter is called the thermodynamic correction factor, , described by The relation between and is then described by equation (14) [175] The collective and transport diffusivities can be calculated from both equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD) simulations. In the latter approach, the chemical potential gradient is the driving force for transport which is imposed on the system in the dual control volume grand canonical molecular dynamics (DCV-GCMD) [199,200].
In EMD, the collective diffusivity can be computed from either of the following equations [149,175] Give the relation of fugacity with chemical potential ∆ ≡ ln ( ), one can rewrite equation (13) in the following form [175,201] ≡ ( ln ln ) (17) here, f represents the fugacity of the bulk fluid in equilibrium with the adsorbed phase and c denotes the concentration of the adsorbed phase.
The thermodynamic correction factor can be calculated from the adsorption isotherm, which itself is obtained from GCMC simulation as explained in Section 6.2.1. Substituting equation (17) in equation (14) will lead to the following relation for the transport diffusivity ( ) [201] ( ) = ( ) ( ln ln ) Diffusion of multi-component systems: To this point, we have discussed methods required for the calculation of different types of diffusion in single-component systems. Diffusion in multicomponent systems is generally an advanced topic with extensive literature available on the fundamentals and practical applications [202]. Here, we mention only essential concepts to illustrate what properties can be obtained from molecular simulations and challenges associated with the incorporation in the process models.
Several equivalently rigorous formulations of multicomponent diffusion exist: e.g., Onsager, Maxwell-Stefan, and the generalized Fick's approach [203,204]. Briefly, for an n-component system, the generalized Fick's law can be formulated as: here, [J] is the column vector of diffusion fluxes of the components in the system and [∇ ] is the column vector of the diffusion gradients in the adsorbed phase. The mutual diffusion matrix, [ ], is given by: with being the Kronecker delta and Ð is the Maxwell-Stefan diffusion coefficients, and [ ] is a matrix of thermodynamic correction coefficients. Equivalently, equation (21) could be formulated using a matrix of Onsager coefficients [ ], which can be shown to be related to [ ] −1 [202,205].
In principle, all properties in equation (20) can be obtained from molecular simulations. Mutual diffusion coefficients and the components of the Onsager matrix can be obtained using expressions, similar to equation (15) for a multicomponent system: whereas elements of [ ] could be obtained from GCMC simulations of multicomponent systems. This immediately points to two challenges. Firstly, construction of the comprehensive data for multicomponent diffusion requires a substantially larger number of simulations, with properties, such as difficult to converge. The complete matrix of thermodynamic correction factors also requires GCMC simulation of multicomponent systems, which may be associated with substantial parameter space (i.e. the variation of the composition of the gas and adsorbed phases). Secondly, the process simulations require a continuous analytical model of the transport and equilibrium properties. Hence, the data obtained from molecular simulations for the properties above would need to be fitted to some simplified models (e.g. the Darken approximation of Maxwell-Stefan coefficients) or be amiable to interpolation within the process model. This will be further complicated, if one wants to incorporate temperature dependence of the diffusion coefficients, since in the micropores it is an activated process.
However, here we do not want to give an impression that the challenges described above present an insurmountable barrier for the implementation of the multiscale workflows for the adsorption problems at hand. In the gas phase, concentration dependent diffusion coefficients would be required for the cases when the number of components is more than two and when the system is expected to significantly deviate from the ideal gas. This is not the case for low pressure binary mixtures of N2 and CO2. Hence, as we will see in the process modelling section, for the diffusion in the gas phase we have a range of classical models that provide concentration independent Fickian diffusion coefficient. In the same section, we will also discuss the fact that in commonly adopted process models for PSA postcombustion carbon capture, the diffusion in micropores of the crystal is not considered at all. The assumption is that for micropores larger than 0.4 nm, the micropores are in instant equilibrium with the gas phase in the macropores of the pellet and we will provide a comment on why it is a reasonable assumption.
Hence, the remaining domain of processes and application where the multicomponent data are indeed required in sufficient detail is associated with kinetic separations, such as for example the separation of oxygen and argon in molecular sieves, or propane-propylene separation using 4A zeolites. Even in the kinetically controlled systems, single component diffusivities coupled with the gradients of the chemical potential will provide a reasonably good model for process simulations in the most cases. Cross-term diffusion coefficients are difficult to parametrise and may not be needed after all. Molecular simulations, however, should provide a possible way to establish when models of intermediate complexity work reliably. We are not, however aware of process modelling studies that incorporated the description of the multicomponent diffusion in its full complexity, although some simple models based on concentration independent single component data have been employed [206,207].

Molecular Dynamics Codes
In this section, we briefly introduce some of the most important open-source molecular dynamics simulation software. There are numerous MD codes developed by various research groups or commercial developers [154,208], some of which are purpose-built software that are developed with particular applications in mind such as those mainly developed for simulation of large biological systems (e.g. NAMD [209], CHARMM [210]). In this section however, we only focus on MD packages that offer many useful features for simulation of fluid transport in nanoporous materials. These softwares are listed in Table 5 and are briefly described here: LAMMPS: the Large-scale Atomic-Molecular Massively Parallel Simulator (LAMMPS) is a highly efficient and scalable classical molecular dynamics simulation code developed by the US Sandia National Laboratories with a focus on materials modelling [211]. It can be used for simulation of solidstate materials (metals, semiconductors), soft matter (biomolecules, polymers), coarse-grained, and mesoscopic systems [211]. LAMMPS can be employed as a parallel particle simulator at the atomic, meso, or continuum scales [211]. LAMMPS is written in C++. Many features of the code support accelerated performance on CPUs, GPUs, Intel Xeon Phis, and OpenMP [211] .

GROMACS:
The Groningen Machine for Chemical Simulations is a MD simulation software primarily designed for simulation of biochemical molecules [212], however due to its computational efficiency it is also highly popular in the domain of materials modelling and simulation of transport processes in porous media. The code is written in C/C++. It was originally developed at the Department of Biophysical Chemistry in University of Groningen. Since 2001, the GROMACS development teams at the Royal Institute of Technology (KTH) and the Uppsala University in Sweden have been responsible for development and maintenance of the software.
DL_POLY: is another classical MD simulation software, which was developed at Daresbury Laboratory in the UK. It is a massively parallel code written in Fortran 90 which is suitable for simulation of macromolecules, polymers, ionic systems, solutions, and transport in porous media [154,213].

Force fields
Atomic force fields for molecular simulations: As discussed earlier in this review, multiscale screening approaches combine various molecular simulation techniques with process modelling and optimization. Molecular simulations are used to compute a series of microscale properties that are required for process modelling including both adsorption data (e.g. equilibrium adsorption isotherms, heat of adsorption, micropore diffusion) and thermal properties of adsorbents (e.g. heat capacity and thermal conductivity). Any inaccuracy in predictions of molecular simulations will manifest itself in the form of errors that will propagate through the process modelling and optimization stage and will lead to inaccurate ranking of porous materials. One key aspect in accurate prediction of molecular simulation data is having access to transferable and accurate atomic force fields. These force fields consist of a set of equations and parameters that describe how molecules interact with each other and with their environment [214] as defined by equation (23) which is a Taylor expansion of classical molecular energy contributions [154,214]: Developing high quality force fields for simulation of different molecular properties across large databases of materials has proven to be a difficult task. In the context of materials screening, the issue with atomic force fields is three folds: (1) availability, (2) consistency and (3) transferability, which will be discussed in this section. For adsorption processes, a chosen set of force fields must at least describe dispersion and electrostatic interactions between all atoms, which are not covalently bonded (assuming that the framework is rigid). Depending on the system, intramolecular interactions of covalently bonded atoms may also be required (e.g. flexible porous frameworks).

Electrostatic interactions:
To compute electrostatic interactions in molecular systems, atomic point charges of all atoms must be known. For adsorption in extended solid frameworks (e.g. porous materials) partial point charges must be assigned to the framework atoms. Although there are multiple computational techniques for partitioning the net electron density (which are normally obtained from quantum mechanical calculations), none of these methods provide an unambiguous definition of the resulting point charges [189]. Unfortunately, an overwhelming majority of material databases constructed so far do not include atomic charges. This can cause a major inconsistency issue in simulation of adsorption processes, which are strongly influenced by electrostatic interactions (e.g. adsorption of polar gases in MOFs and zeolites). Effect of framework atomic charges on gas adsorption have been already demonstrated in several studies [215,216]. It is shown that application of different charge calculation schemes in molecular simulation can easily become a source of inconsistency in comparing various screening studies [215,216]. In this context, Li et al. [41] provide a good example for the use of two different charge calculation techniques for screening of MOFs for CO2 capture in the presence of water. They use the extended charge equilibration (EQeq) method as well as the repeating electrostatic potential extracted atomic (REPEAT) method to compute atomic partial charges of MOFs for adsorption of CO2/H2O and CO2/H2O/N2 mixtures in a large group of MOFs, which are selected from the CoRE-MOF database. It is demonstrated that electrostatic interactions play a dominant role in determining the water adsorption behaviour in MOFs and that a judicious decision must be made in choosing the method for calculating framework atomic charges [41]. In a more recent study, Altintas and Keskin [217] discussed the role of partial charge assignment methods in highthroughput screening of MOFs for CO2/CH4 separation. They employ a quantum-based density-derived electrostatic and chemical charge method (DDEC) and an approximate charge equilibration method (Qeq) to predict adsorption of CO2/CH4 mixtures in 1500 MOFs. The authors show that gas uptake, working capacity and mixture selectivity of MOFs vary considerably depending on the charge assignment methods used in molecular simulations. They also report that the ranking of the bestperforming MOFs are also different depending on the method used for charge calculations [217].

Dispersion interactions:
Another important concern regarding the use of atomic force fields is availability and transferability of model parameters for calculation of dispersion interactions (e.g. Lennard-Jones parameters). Currently and for practical reasons, high-throughput screening of large materials databases heavily rely on the use of generic force fields [39-41, 48, 55, 109, 218, 219]. The most commonly used generic force fields include DREIDING [220], UFF [221] and OPLS-AA [222]. Despite their wide spread, generic force fields fail to accurately reproduce experimental adsorption data in many cases [223][224][225] particularly for gas adsorption in MOFs with coordinatively unsaturated metal sites [226][227][228]. Even for the systems where generic force fields are deemed suitable, prediction of experimental adsorption isotherms is rather qualitative in which simulated isotherms only capture the general shape of their experimental counterpart [229][230][231][232]. Therefore, the use of generic force fields for screening of large and diverse databases of porous materials should be approached with caution [229,233,234]. In fact, many researchers have already started to develop specialized force fields for challenging systems such as those involving adsorption of water [197,225] or MOFs containing open metal sites [227,[235][236][237][238].
For adsorption, normally a set of force field is considered adequate if it can reasonably reproduce experimental adsorption isotherms. Nevertheless, one major issue in quantifying the adequacy of atomic force fields is associated with what is perceived as a "good agreement" between experimental and simulated adsorption isotherms. In many studies reviewed above, simulated adsorption isotherms obtained based on the use of generic force fields only qualitatively capture the overall shape of the experimental isotherms. Although this may seem sufficient from a molecular simulation point of view, it is apparently inadequate for accurate prediction of adsorbents performance in the process level. This is particularly shown for PSA/VSA processes, which are reviewed in this article [50,51,83]. These studies have demonstrated that separation performance of adsorbent materials (i.e. their energyproductivity Pareto fronts) is sensitive to key characteristics of mixture adsorption isotherms such as Henry's constant, non-linearity and selectivity [50,51,83].
Another important concern regarding the quality of currently available force fields is associated with the role of lighter components (e.g. nitrogen) in process simulation. To date, significant efforts have been made to develop more accurate force fields for adsorption of target gases in important separation processes (e.g. CO2 in carbon capture, CH4 in natural gas purification), while the importance of the other components in the mixture has been somewhat overlooked. In the context of postcombustion carbon capture, almost the entire effort of molecular simulation community has been dedicated to the development of more accurate force fields for CO2. However, from a process simulation perspective, adsorption of nitrogen plays an important role in separation performance of adsorbents in PSA/VSA processes [6,50,51,55,79,84,239]. Therefore, it is equally important to develop accurate and transferable force fields for simulation of nitrogen adsorption in the diverse group of porous materials, which are currently available from various databases. A quick review in the literature shows that specialized QM-derived force fields for nitrogen adsorption are scarce [240][241][242][243], in spite of the fact that performance of generic force fields for prediction of nitrogen adsorption in many materials is not satisfactory and does not meet the level of accuracy required for process simulation. Some examples of this include adsorption of nitrogen in STT, CHA all silica zeolites [244], FAU and MFI type zeolites with different Si/Al ratio [51,244], Mg-MOF74 [240] and Ni-MOF74 [83] , ZIF-68 [245], Zn-MOF and Cu-BTC [246]. It is therefore important to develop a cross-disciplinary awareness about the implications of the accuracy of atomic force fields which are developed on a molecular level, but their adequacy will have profound impact on the process level predictions in the materials screening studies.

Force fields for prediction of atomic vibrations:
Molecular simulations typically assume adsorbent materials to have rigid atomistic frameworks. Recently, novel porous materials have been discovered that exhibit structural flexibility [247,248]. Development of force fields that can correctly capture this behaviour is an ongoing area of research. This is particularly important for the studies of MOFs, as all MOFs exhibit some form of structural flexibility [157,214,247] ranging from lattice vibrations at equilibrium to large-scale structural transformations upon external stimuli [214], such as temperature [249], guest adsorption [250] and electric field [251]. Among different types of structural flexibilities, structural vibrations and phonon properties of the lattice determine specific heat capacity of porous solids [190,252,253], whose relevance for process simulation has been recently explored [83].
As elucidated by Kapil et al., thermal properties of the lattice can be described by a quantum harmonic treatment [253]. However, the heat capacity of loaded porous frameworks requires a combination of quantum and anharmonic treatment [253]. Analysis of phonon properties for estimation of thermal properties of materials requires costly quantum mechanical calculations [254] which are not affordable for routine screening of large numbers of porous materials. To address this limitation, development of purpose-built and computationally affordable force fields has been recently undertaken by several groups [254][255][256][257], nevertheless further developments for improved accuracy and transferability of these force fields are required [258].

Process Modelling and Optimization
The main objective of this section is to give an accessible guide on the PSA/VSA process modelling from fundamentals to practical implementation. We begin with the basics of the mass, energy and momentum balances in the adsorption column packed with pellets of adsorbent material (Section 6.3.1). We will introduce the hierarchy of models, differing in the levels of details in their description and in the assumptions involved. We will briefly review the commonly involved methods in the solution of the introduced balance equations, under the appropriate boundary conditions.
Setting up a process model requires a number of parameters and properties. For a non-practitioner, it can be overwhelming to see the process model in its full complexity, and hence in the next Section 6.3.2 we tasked ourselves with explaining what parameters are required and how their values can be obtained.
A pressure swing adsorption process involves several columns, each of them going through a cyclic sequence of steps. In the next Section 6.3.3, we will use a simple 4-step cycle to introduce the PSA process and the key concepts associated with its cycle, such as cyclic steady state (CSS), performance of the cycle in terms of purity, recovery, productivity and energy consumption. Furthermore, using this example of the 4-step process, we will briefly explore the concentration profiles during different steps at CSS and how to interpret them.
A specific cycle configuration may not operate at the optimal conditions. Hence, a significant part of the process modelling research is focused the cycle optimization. In the next section 6.3.4, we introduce currently used optimization methods, such as genetic algorithms, and essential concepts associated with the process optimization.
Initially we had focused on the example of a 4-step cycle, more complex cycles involving larger number of steps and columns as well as multiple stages have also been explored in the context of postcombustion carbon capture. We believe it is useful to review these studies and the drive behind the exploration of the more complex processes. This is addressed in Section 6.3.5.
As has been already discussed in the section on the process metrics, in general process simulations are time consuming. This prompted a significant research effort into more efficient alternatives for process performance evaluation that work in tandem with detailed process simulations. These developments are reviewed in Section 6.3.7.
Finally, following the spirit of the review, we conclude the section on process modelling with a brief overview of the available codes for this type of modelling, their capabilities and access (Section 6.3.7). An adsorption column is the basic unit of the adsorption process. In this section, we provide a brief summary of the mass, energy and moment balances around this unit, which are either solved numerically in the process simulations or serve as starting points for simplified analytical models. For a more comprehensive analysis we refer the reader to the seminal books by Ruthven et al. [259,260] on fundamentals of adsorption and PSA processes.

Fundamentals
Consider a schematic of a packed column in Figure 9 . The column has length of , is used as the position within the column in the axial direction and the feed is introduced to the column from the bottom at = 0. The column is packed with pellets of adsorbent material. The pellet consists of microporous crystallites, which are held together by an inert binder. Thus, the pellet has intercrystalline macropores and intra-crystalline micropores. In the description of the various transport processes, we adopt the following convention: macropore refers to the pore space between the crystallites and micropore refers to the pores inside the crystallites. On the right of Figure 9 , we show an idealized spherical pellet of radius and volume . In the model, we can also assume that crystallites are spherical particle of radius . The pellet volume consists of the macropore volume and crystal volume , which in turn consists of the micropore volume and the skeletal volume : The bulk density , pellet density , crystal density and skeletal density are defined as follows: Here, is the volume of the gas phase in the column and is the total mass of the adsorbent pellets. This mass includes both the mass of the adsorbent crystals as well as the mass of the binder. Thus, it is assumed that the binder volume is part of the skeletal volume of the pellet. Therefore, the saturation capacity of the adsorbent has to be corrected for the mass of the binder if the adsorption isotherms were measured for the non-pelletized adsorbent crystals. The bed void fraction , pellet void fraction and crystal void fraction are defined as follows: A common starting point for many process modelling approaches is the material balance in the column based on the axial dispersed plug flow model with an adsorption source term (although more complex and complete formulations are also possible, i.e. including radial dispersion term, etc): Here, is the gas phase concentration of component i, is the macropore concentration of component i in the adsorbent pellet and is the sorbate concentration of component i in the adsorbent, is the interstitial velocity, and is dispersive flux of component i. In this equation, the first and the second terms are the accumulation terms in the gas phase and in the pellets, respectively. The amount adsorbed in the pellet can be seen as the composite of the amount as gas in the macropores of the pellet and the absolute amount adsorbed in the micropores of the adsorbent material, (1 − ) : In the column mass balance, the average amount adsorbed in the pellet is needed: and, similarly the average adsorbed in the crystallite can be defined: The third term in equation (33) describes the convective flow of the gas across the bed and the final term describes the dispersion process relative to the bulk flow. The dispersive flux is given by: where, is the axial dispersion coefficient, is the total gas concentration, and is the mole fraction of component i. For the axial dispersion coefficient ( ) correlations are available [260], such as: here, is molecular diffusivity which is defined by equation (42), and 0 is the average superficial fluid velocity through the packed bed.
Equation (33) provides the overall mass-balance in the column, it does not describe the actual process of diffusion into the pellets. For this, a separate set of material balance equations can be formulated around the pellet. In the most general case, the model will contain terms associated with the external film resistance at the pellet surface, macropore diffusion from the bulk gas phase into the pellet, barrier and film resistance at the adsorbent crystals boundary and micropore diffusion in the adsorbent crystals. A schematic of an adsorbent pellet with relevant properties is shown on the right of Figure 9. Let us consider these processes in more detail.
First, let us focus on the overall material balance for the pellet. The concentration in the macropores is governed by the following mass-balance equation, based on the second Fick's law formulated for the spherical pellet geometry: here, , is the effective macropore diffusion coefficient. The first term of equation (39) represents the accumulation in the macropores, the second term describes the accumulation in the micropores and the last term describes diffusive mass-transport due to the concentration gradients inside the pellet (the second Fick's law). Here cross interactions are neglected and this is typically sufficiently accurate, especially considering that in most optimization studies an equivalent LDF lumped coefficient is used to describe the mass transport as it captures the essential features at a limited computational cost.
At the surface of the pellet, diffusion from the bulk gas phase into the pellet can be described via masstransfer process across the film at the surface: where , is the external fluid film mass transfer coefficient. This equation sets the boundary condition at , whereas at = 0 the boundary condition is | =0 = 0 which is required due to the assumption of spherical symmetry.
The effective diffusion coefficient reflects various mass-transfer mechanisms into the pellet and is obtained by combining the molecular diffusion , Knudsen diffusion , surface diffusion and the viscous diffusion coefficients : where, the individual diffusion coefficients are estimated using well-known expressions: here, is the molecular weight in g mol -1 , 12 is the collision diameter from the Lennard-Jones potential in Å, Ω 12 is a function depending on the Lennard-Jones force constant and temperature, is the mean macropore radius in m, is the Henry's constant of adsorption, and is the diffusional activation energy. These expressions along with the theories behind them and the values of the parameters are discussed in the classical textbooks on transport phenomena [259,261]. We further note, that typically in the process models the values are obtained at some fixed, representative conditions, while in reality the conditions change dynamically in the actual process, and hence, these properties would also vary in time in a more accurate model.
Similarly, for the diffusive process in the micropores inside the crystallites, modelled as spherical particles of size , we can formulate a similar general mass-balance equation, based on the second Fick's law of diffusion: here, is the effective diffusion coefficient in micropores and other terms are described as before. Similar to the processes at the pellet surface, the diffusion into the crystallite particle from the surface can be described using transfer resistances across the surface: Here, * is the adsorbed concentration of component i in equilibrium and is the concentration of component i at the crystal boundary. In equation (47)  A similar hierarchy of equations can be formulated for the energy balance in the column. In the most general non-isothermal case, the following equations govern heat-transfer processes: here, ̌ is the internal energy in the fluid phase per unit volume, ̌ is the internal energy in the pellet per unit volume, ̌ is the enthalpy in the fluid phase per unit volume, is the thermal diffusive flux, ̃i s the partial molar enthalpy of component i in the fluid phase. , are temperatures of the fluid and the wall, respectively, while ℎ is the heat transfer coefficient between the wall and the surroundings. and are the surface area and the volume of the column, respectively. The first two terms in Eq. 39 are accumulation terms for the gas phase and the solid phase respectively; the third term is associated with the convective flux of fluid stream, with enthalpy ̌. The next two terms are axial dispersion terms. The first one, , describes thermal flux due to the temperature gradients along the z axis, whereas the second term is associated with the diffusive fluxes due to the concentration gradients along z axis (and hence enthalpy fluxes coupled with them). The last term on the left in equation (48) describes heat transfer across from fluid to the wall of the column.
The heat transfer across the wall of the column can be described as: where ∞ is the temperature of surroundings. The column wall is defined by the column wall density and specific heat capacity, ̂, . The ratio of the logarithmic mean surface area to volume of the column wall are given by: where and are the radius of the column and the thickness of the wall, respectively.
At the level of the pellet, uniform temperature profile is typically assumed across the pellet (no temperature gradients) and this has been shown to be consistent with the experimental observations [262].
The overall energy balance for the pellet can be then formulated as: Finally, the thermal axial dispersion flux is given by: Here, the axial thermal conductivity in the fluid and pellet are given by and , respectively. There is also an alternative way to formulate the energy balance for which we refer the reader to the original articleby Zhao et al. [263].
The momentum balance is described by the Ergun pressure drop equation: where, is the fluid viscosity and is the fluid density. Figure 10. Hierarchy of the models available for the mass and energy balance in the adsorption column (not exhaustive list). Squares shaded green reflect the combination of the models employed in the studies by Farmahini et al. [51,83] and also commonly adopted by other practitioners in the field.
The equations above provide a complete and general description of the mass and energy balances in the column. These equations serve as a starting point for more simplified models. Indeed, Figure 10 illustrates the hierarchy of the models with each model based on its own set of assumptions and resulting simplifications of the governing equations. Reading this diagram from left to right, the system can be considered as isothermal (hence no energy balance equations are required) or non-isothermal. Then, within each branch, we can either include or ignore the pressure drop across the system. For each branch we can further consider whether we include film resistance at the surface of the pellet or not and so on. This hierarchy demonstrates that we can construct order of 10 2 models depending on the combination of the assumptions we use. The boxes shaded green in Figure 10 represent the choice of the assumptions adopted in the studies of Farmahini et al. [51,83] , as well as in many other previous studies [78,264,265]. In this case, the following assumptions are considered: 1. The system is modelled as non-isothermal with heat transfer allowed between the packed bed and its wall, but pellets and gas phase are kept at the same temperature. = 2. Pressure drop is considered across the bed. The pressure drop is modelled using the Ergun equation (53).
3. No external film resistance is considered. In this case, Eq. 35 vanishes and the following condition applies: ( , , ) = ( , ) 4. The macropore resistance is modelled using the Linear Driving Force (LDF) approximation. Effectively, all the resistances to diffusion are lumped into a single effective parameter, while the driving force of the process is simply the difference between the concentration of species i in the gas phase and in the macropores, . As a result, Eq. 34 can be replaced with a simplified model: here, is the LDF coefficient for the pellet, which can be calculated using effective macropore diffusivity with the Glueckauf approximation, which is equivalent to assuming a parabolic concentration profile [266]: 5. Micropore equilibrium is assumed. This assumption implies that the crystallites are in instant equilibrium with the gas phase in the macropores of the pellet. This would be the case when the overall mass-transfer into the pellets is controlled by macropores, and not micropores. Although this seems counter intuitive, the validity of this assumption for materials with pore sizes > 4 Å has been discussed on several occasions [259,267]. Briefly, the reason for this is as follows. Although the diffusion coefficient is several orders of magnitude lower in the micropores of the crystalline material compared to the diffusion in the gas phase in the macropores, the diffusion in macropores happens on the scale of the pellet, and hence involves longer characteristic time scale. To illustrate this point, let us return to the equation (39) describing the mass balance around the pellet. If we make an assumption that the isotherm is linear (̅ = , where , is the Henry's constant for component i), equation (39) can be rearranged as: The two diffusional time constant which should be compared to each other are then the macropore diffusion constant 2 / , and the micropore diffusion time constant, 2 / . While it is obvious that is always smaller than , , what is important in determining the controlling mass transfer mechanism is the comparison of the molar fluxes, and often macropore diffusion is the controlling mechanism due to the combined effect of small crystals in relatively large beads and the large value of the effective bead Henry law constant. This assumption is equivalent to the following condition: If instead of the condition above, one may wish to include a more detailed model of micropore diffusion using the LDF approximation, then the following simplification can be employed to describe transport into the crystallites: where is the combined LDF coefficient for the adsorbent crystal, and is the LDF coefficient for the micropore; , is the mass-transfer coefficient across the external fluid film, and , is the masstransfer coefficient at the crystal boundary defined earlier. For further details on grouping transport resistances into LDF-equivalent coefficients we refer the reader to the existing reference books on fundamentals of adsorption and PSA processes [259,260]. The LDF coefficient can be calculated from the effective micropore diffusivity by: Regardless of the details of the model, the combined mass and energy balance equations form a system of Differential Algebraic Equations (DAEs). These equations are usually discretized in the spatial domain by an appropriate numerical method such as finite difference, finite element, orthogonal collocation, or finite volume method. This produces a system of Ordinary Differential Equations (ODEs) which can be solved using a number of approaches: either internal functions within the available simulation packages, such as MATLAB, or existing numerical solution libraries (e.g. SUNDIALS [268] ).
Here, it is also useful to reflect on the simplified LDF-based models versus detailed diffusion equation models. The motivation to develop simplified LDF models is driven primarily by the numerical efficiency. Indeed the simplified LDF model with 30 axial volumes and 2 components corresponds to 120 DAEs, while the same with the diffusion equation would be approximately 600. Including also diffusion in the micropores would be ca. 3000 equations. The computational costs would be at least proportional to the total number of equations N if the code is well written and N 2 for a not so well written code.

Complete Hierarchy of Data Required for Multiscale Process Simulation
One of the primary aspirations of this review is to provide a useful guide on PSA/VSA process models for non-practitioners. Reading a standard research paper on process modelling of adsorption processes can often be overwhelming because of the number of parameters and properties one needs to specify, with their sources not necessarily being obvious. Here, we also emphasize that even after reading our review we do not expect a novice in process modelling to be able to setup their own simulations. However, we hope they will be able to understand the requirements for these simulations and critically appraise their source, reliability and relative importance. Broadly, we can split the data required for setting PSA/VSA process simulations into the following categories. Column properties describe the geometric dimension of the column, its length, diameter, and thickness of the walls. These properties either are taken to reflect the actual experimental unit, or given some specific, physically meaningful values. For example, certain parameters of the column have been used in several studies and they have now become commonly employed parameters for several groups to ensure consistent comparison of the process modelling results in different studies [50,51,55,77,[81][82][83][84]. The balance equations described in section 6.3.1, also imply that to solve these equations we need values of the properties associated with the thermophysical characteristics of the material of the column and how it interacts with its environment (e.g. heat capacity, heat transfer coefficient, etc).
In the next category, we have all properties associated with the pellet: pellet size, pellet porosity, and pellet tortuosity. In the same category, we also include properties associated with the transport in macropores of the pellets, such as different contributions to the overall macropore diffusivity (e.g. molecular diffusion, Knudsen diffusion, etc). Figure 10 is crystallites, and hence the next category of properties is associated with the properties of the adsorbent material crystals: crystal density, crystal thermal and transport properties, etc. In principle, the pellet is made out of crystallites and binder and properties of the pellet, such as the specific heat capacity or thermal conductivity, are a composite property of the two materials, binder and crystallites. However, the common convention is to assume these properties of the binder to be equivalent to the properties of the adsorbent crystals.

Further down in the hierarchy of scales shown in
Generally, equilibrium adsorption data should also belong to the category of the crystal properties. However, this requires a special consideration. Adsorption data, both in experiments and in simulations, are typically obtained as single component adsorption isotherms comprised of discrete data points. However, process simulations require an analytical expression describing adsorption equilibria in order to be able to solve the mass balance equations described in Section 6.3.1. Moreover, the accuracy of process modelling also depends on how well the supplied models describe the multicomponent equilibria -hence accurate interpolation of single component isotherms may not be sufficient for the correct behaviour of the model in the actual process simulations. A common approach is to use the dual-site Langmuir (DSL) adsorption model to obtain an analytical description of adsorption isotherms. For a single component system, the DSL isotherm for species i is defined by: here, , is saturation capacity of site j with respect to species i, bj,i is affinity of each site described by the van't Hoff equation: , = , exp ( − , ). In the van't Hoff equation, , is the heat of adsorption at adsorption site j and , is a pre-exponential factor.
As seen here, there are six parameters ( 1, , 2, , 1, , 2, ,  1, and  2, ) for each gas component i which can be obtained. The thermodynamic consistency requires that the saturation capacity of each site is the same for all adsorbing species (for example, for the binary CO2/N2 adsorption this implies 1, 2 = 1, 2 and 2, 2 = 2, 2 ), unless adsorbing molecules differ significantly in size. In an early study, Myers showed that these conditions are essential for the accuracy of the multicomponent DSL model [269]. This poses additional constraints on the fitting of equation (65) to the reference adsorption data using non-linear least-square regression. Adsorption of species A from a binary gas mixture of A and B at fixed temperature is described by the extended version of the Dual-Site Langmuir model (extended DSL) which is given by: where yA and yB are mole fractions of components A and B in the gas phase. To obtain physically meaningful parameters for the DSL model, normally the fitting algorithm is guided through a set of mathematical constraints, which also help the algorithm to converge [51]. The quality of the DSL model is ultimately tested by its ability to predict binary adsorption equilibria. This data may not be readily available from experiments, however in molecular simulations it is relatively easy to implement and carry out these tests. In our previous publications, we explored systematic ways to obtain parameters of the DSL model and we refer the reader to the original publication [51].
It should be noted that for many new materials such as phase-change adsorbents it is not easy to propose a suitable functional form that can properly describe equilibrium adsorption data [270,271]. The alternative approach here is to describe the equilibrium relationship between adsorbed phase and fluid phase as a set of discrete points. Haghpanah et al. [272] have proposed a method to obtain discrete equilibrium data from single-component breakthrough experiments and include it into computer simulations so that a continuous functional form is no longer required. In this method, adsorbed phase concentration (q) is defined for a set of discrete values of the fluid phase concentration (c) within the range of the feed concentration. The adsorbed phase concentration of any point between two adjacent discrete points is then calculated by interpolation [272]. Finally, the adsorbed phase concentrations is computed by solving an optimization problem where for a given set of discrete fluid phase concentrations, the experimental breakthrough profile is reproduced [272]. The above computational technique has been further developed by other research groups [273,274]. For example, Rajendran et al. [274] have extended this method by incorporating discrete singlecomponent equilibrium data into the Ideal Adsorbed Solution Theory (IAST) [275] in order to describe binary equilibrium data.
The final category of parameters that are required for process modelling include properties of the feed such as its temperatureand composition which are typically specified by the design problem at hand (e.g. post-combustion carbon capture). Table 6 summarizes the full set of properties needed to set up a PSA/VSA process simulation along with their sources according to the categories provided above. From the table above it is clear that setting up a model requires a combination of properties that can be measured experimentally (e.g. adsorption isotherms, properties of the pellet), or for which wellestablished thermophysical models exist (e.g. molecular diffusivity, Knudsen diffusivity). Some other properties have well-known literature values (e.g. heat conductivity of steel). In general, the large number of parameters required to set up the model in combination with the large number of potential models (hierarchies used as described in Figure 10) often make comparison and reproduction of data between various research groups a challenging task and we can only advocate detailed disclosure of the sources, parameters and algorithms used for every simulation.
A separate challenge is the implementation of the complete in silico workflows. As can be seen from Table 6, only a limited set of properties can be obtained from molecular simulations (e.g. equilibrium adsorption data, micropore diffusivity, heat capacity and thermal conductivity of adsorbent crystals). For other properties, particularly those pertaining the morphology of the pellets, we can either adopt some conventional estimates based on what is known from previous experimental measurements or we can use these parameters as optimization variables within a specified range of known values. The former approach is however prone to inaccuracy and inconsistency, considering pellet morphology is not standardized and various manufacturers produce adsorbent materials with different characteristics (e.g. different size and porosity, various types of binder). Optimization of these parameters however have proved to be a more promising approach in some cases. In a recent study, Farmahini et al. [83] have demonstrated that size and porosity of pellets can be used as decision variables during process optimization not only to achieve maximum theoretical performance of adsorbent materials, but also for consistent comparison of different materials screening studies [83].
To fully understand the impact of these two approaches, we advocate for sensitivity and error propagation analyses of the multiscale materials screening workflows for the parameters that cannot be calculated from molecular simulations or any other theoretical methods. The results of such analyses will show whether the use of estimated reference values for these properties has a significant impact on the overall predictions of the multiscale workflows.

PSA/VSA Process and Cycle Configuration
In Section 3, we briefly introduced the PSA/VSA process. In the previous sections, we also covered the mass, energy and momentum balance equations, governing the behaviour of the adsorption column and data needed to set up the process model. Here, we consider in more detail a particular 4-step VSA cycle and essential elements of cycle configuration. For the sake of concreteness and consistency we continue with the same case study of the post-combustion feed, comprised of carbon dioxide (15%) and nitrogen (85%).  [84]. Figure 11 shows a 4-step cycle which first appeared in the work of Ko et al. [278], who referred to this as the fractionated vacuum swing adsorption cycle. The first step of the process is the adsorption step. The feed is introduced to the column at the pressure close to atmospheric. This is followed by a cocurrent blowdown step: the column is closed at the feed end, and the pressure is reduced to remove excess nitrogen present in the column in order to increase the purity of the product. Next is the counter-current evacuation step, where the pressure is reduced further, causing desorption of carbon dioxide. The product of this step is a carbon dioxide-rich stream. Finally, this step must be followed by bringing the pressure of the column back to the adsorption pressure, which is done in the repressurization step. In principle, repressurization can be done using the feed stream. However, previous studies demonstrated that counter-current repressurization with the light product stream, as schematically depicted in Figure 11 leads to a much better process performance [80,279]. This effect stems from counter-current repressurization helping to concentrate carbon dioxide closer to the feed end of the column as it will increase purity and recovery of carbon dioxide during the evacuation step later in the sequence. As can be seen from Figure 11, each step in this process is therefore associated with a particular pressure profile and duration. These parameters, namely time of the adsorption step , time of the desorption step , and time of the evacuation step , the blowdown pressure , and the evacuation pressure , along with feed pressure , and feed flow rate are called cycle variables for this particular process and their specific values define the cycle configuration. The cycle variables are typically constrained by a number of considerations. For example, cannot be set too high otherwise the compression of the dilute gas makes the process not viable.
is another important example. In practical systems 0.2 -0.3 bar would be a reasonable value for this parameter, but often much lower values are used in process simulations in order to achieve the required purity/recovery targets. As a case study, let us consider simulation of a single cycle configuration using the values of the cycle variables specified in Table 7. The values of other properties required to set up this simulation have been specified in the previous section. Adsorption processes operate at the Cyclic Steady State (CSS), and equations described in 6.3.1 can be solved iteratively to arrive to the CSS. Alternatively, time can be discretised and the CSS in this case is calculated directly, but this approach resulting in a large set of nonlinear equations is not necessarily faster [280]. Although the actual industrial process features several adsorption units in a different stage of the cycle at any given moment, as they all go through the same steps, it is possible to consider modelling of this process with only one unit. This so-called unibed approach has been originally described by Kumar et al. [281]. This in general allows to study multicolumn process at a similar numerical cost as a simple Skarstrom cycle. The solution procedure starts with some initial conditions and solution of the balance equations in the adsorption step. This produces concentration profiles for each component of the system in the adsorbed phase and in the gas phase. These concentration profiles and the composition of the product stream serve as the initial conditions for the next step in the adsorption cycle (in this case, the blowdown step), and so on. The iterative process continues until the numerical CSS is reached: this happens when the state variables start to depend only on the spatial position in the system and the time relative to the start of the cycle. One can employ several mathematical criteria to establish whether the solution has reached the CSS [282,283]. We note, however, that this is not a simple problem especially for non-isothermal systems and (or) with one very strongly adsorbed component (for example water): in this case convergence may requires thousands of cycles [283]. Complexity of the PSA/VSA processes is very well illustrated by the concentration, temperature and pressure profiles that are calculated for each process cycle. Correct interpretation of these data is vital for the analysis of the performance and efficiency of the PSA/VSA processes.  The CSS implies that these profiles do not change anymore (within the numerical convergence criteria) as we continue with the numerical iterations and they will remain looking like this at the end of their respective steps.
Let us focus on these profiles in a step by step fashion. The LPP step prepares the bed for the next adsorption step and the LPP profile reflects the state of the column before the adsorption step is started. Compared to the process not featuring LPP, the main advantage of LPP is that this step pushes carbon dioxide towards the column feed and this generally leads to the improved recovery. In the gas phase the concentration of carbon dioxide is very low. In the adsorbed phase, the concentration of carbon dioxide is also low, however some carbon dioxide remains in the adsorbed phase close to the feed end of the column. At the end of the adsorption step, the profile in the adsorbed phase reflects the higher amount of carbon dioxide now present in the solid. It starts with saturation values at the feed end slowly diminishing towards the light product end of the column. This reduction in saturation value is due to a non-uniform temperature distribution along the column: at the adsorption front the heat of adsorption increases the temperature which in turn reduces the saturation value; behind the adsorption front the temperature reduces slowly and, in turn, the saturation value increases towards the feed end. In the gas phase, the concentration of carbon dioxide is low at the end of the adsorption step. The main purpose of the blowdown step is to remove remaining nitrogen in the gas phase. At the end of the blowdown, some of carbon dioxide is released from the porous material and it is concentrated at the feed end of the column in the gas phase. The available carbon dioxide at the end of the blowdown step will contribute to the heavy product, i.e. CO2-rich stream, during the countercurrent evacuation step. At the end of this step, the gas phase consists almost of pure CO2, while in the adsorbent phase the concentration of CO2 is lowered. It is important to note from the profiles discussed that the porous material is never fully regenerated -the amount of carbon dioxide it captures is represented by the difference between blue (adsorption) and the green (evacuation) lines, indicating that the working capacity of the material is only a fraction of the absolute capacity. From the same graph it is also clear that the adsorption step is stopped before the complete breakthrough occurs and the portion of the bed (between 0.75 and 1.00 in the dimensionless coordinates along the bed length) is never used.
To quantify performance of PSA/VSA processes, the following properties are normally evaluated: 1) Purity, 2 : this property characterizes the composition of the final product. It is the ratio of the number of moles of carbon dioxide evacuated to the total number of moles of gas mixture evacuated during a single cycle:

=
Moles of CO 2 recovered in evacuation Total moles out in evacuation 2) Recovery, 2 : this property describes the amount of carbon dioxide recovered as part of the product stream compared to what is originally fed into the column.

Moles of CO 2 recovered in evacuation
Total moles CO 2 in the feed The other two properties include energy penalty and productivity of the process, which have been already defined in Table 1, however we will explain them here again: : 3) Energy penalty: it is defined as the total amount of energy used for separation of one mole of CO2 from the feed.
here, it is also instructive to reflect on the nature of the energy used in the process. In the PSA/VSA cycle this work is associated either with compression or pulling vacuum. In the 4-step process considered here, the most significant energy penalty comes from pulling vacuum during the evacuation step, however it may shift to other steps in more complex processes [284].
The complexity of this picture, its dynamic nature and the fact that it depends on a number of parameters, including the configuration variables of the cycle, explains why it is difficult to find some simplified metric which would comprehensively capture efficiency of PSA/VSA separation processes.

Process Performance and Optimization
In the previous section, we considered a single cycle configuration with specific values ascribed , , , , , and flow rate of the feed, F. However, in reality, the resulting process may or may not be able to meet the design objectives to recover more than 90% of CO2 with at least 95% purity. It also may not operate optimally; hence incurring additional energy penalties. The objective of the optimization process is to adjust the values of the cycle parameters in such a way that the process can meet its design constraints, while operating at the highest possible productivity and minimum energy penalty. In the optimization language, the cycle parameters described above become decision variables, while mathematically the optimization problem can be formulated as follows: The optimization conditions above form an optimisation problem with two objective functions and two constraints. Here, it is important to realize that the two optimization targets, minimal energy penalty and high productivity, are in competition with each other. Indeed, higher productivity may be achieved using higher flow rates given the same amount of the active adsorbent material in the column. However, this approach may require faster cycles and lower evacuation pressures, which will lead to higher energy penalties. In contrast, lower energy penalty can be achieved with more moderate vacuum during the evacuation step but it will be achieved at a cost of processing lower flow rates in the system or having to resort to longer individual steps, leading to lower productivity. As a result, the actual solution to the optimization problem is not a single set of values of the cycle parameters, but multiple combinations of these parameters, each of them associated with a particular combination of purity, recovery, energy penalty and productivity values.
From the mathematical perspective, the problem above corresponds to a multi-objective optimization. In general, this is a challenging problem as the search for the solutions takes place in a multidimensional space of the decision variables, which can form clusters of feasible solutions, separated by non-feasible regions. The study of Fiandaca et al. [285], showed that the objective function is non-smooth and non-convex, and also that the design space is non-convex. Several approaches have been proposed to deal with this problem over the years, with Ref [285] briefly reviewing available approaches up to 2009. However, in recent years the conventional practice became to invoke the evolutionary Genetic Algorithms (GA), because of their ability to achieve global convergence, and a large number of tools available to implement them. In particular a set of methods associated with the second and third generation of non-dominated sorting genetic algorithm (NSGA-II,III) has been a popular choice. It has been implemented in many commercial packages such as MATLAB and also available as a set of free libraries [286][287][288].
The initial step in the optimization problem is to identify a range of values within which each decision variable can change. A number of initial operating conditions (so called, population in GA terms) is selected from this range (either randomly or using more sophisticated approaches such as Latin Hypercube Sampling). For each combination of the decision variables, the PSA process is simulated as described in Section 6.3.3. Promising candidates are identified, and their features are combined (using mutations and crossover moves) to give a new generation of operating conditions. As the optimization process evolves from generation to generation, the cloud of points representing the cycle configurations on the Energy Penalty-Productivity progresses towards higher values of productivity and lower values of energy (subject to purity and recovery constraints) until this process effectively stops (further progress of the cloud is not visible within the convergence criteria). At this point, the optimization simulation has converged to its final set of solutions.
This process can be illustrated with two useful graphs commonly employed in the process simulation and optimization studies. The first plot has purity and recovery as X and Y axes. It identifies the proportion of cycle configurations that are able to meet the 95%-90% constraints for purity and recovery. The second plot shows the evolution of the cycles in Energy Penalty -Productivity coordinates. Figure 13 illustrates typical examples of these graphs. Figure 13. Process performance characterized in terms of Purity-Recovery coordinates (constraints, top graph) and Energy Penalty -Productivity coordinates (Pareto front, bottom graphs); graphs on the left is a schematic for the illustration; whereas graphs on the right correspond to a case studied in our recent publication [83].
The front edge of the clouds shown in the above figure are called the Pareto fronts. These are the set of cycle configurations that combine the highest purity-recovery and energy-productivity for a given process configuration subject to its pre-defined process constraints.
As already mentioned, this implies that for each material there is a number of possible operating conditions to choose from (points on the Pareto front). High productivity processes will incur higher energy cost, but lower footprint and capital cost of the plant. Low energy processes, on the other, hand will benefit from lower energy penalties, however, may incur larger capital costs due to larger required footprint of the plant. these values correspond to: to the lowest energy penalty on the Pareto front, or the highest productivity.

More Complex Process Configurations and Recent Studies
The main objective of this section is to introduce the reader to more complex PSA/VSA process configurations and review recent studies on application of process modelling to assess the viability of PSA/VSA technologies for carbon capture.
In the previous section, we used the 4-step VSA-LPP process to introduce several essential concepts and fundamentals of the PSA/VSA process and optimization. One of the issues associated with this specific process is that it can meet the required purity/recovery constraints only by going to very low evacuation pressures (e.g. 0.01 bar). Although from the Pareto front analysis this process is very competitive compared to other alternatives, in practice it is not viable, as the standard industrial pumps do not typically go below the range 0.13-0.2 bar [289]. This necessitates a search for more complex process configurations. One option is to consider a two-stage process [267]. In this case, the first stage focuses on maximizing the recovery, while the second smaller polishing unit would aim to achieve the required purity. Indeed, Abanades et al., summarized recent studies of 2-stage PSA processes [267]. According to their summary, it is clear that most process simulations arrive at VSA configurations that require approximately 0.5-0.75 MJ/kg CO2, and that they can operate at a evacuation pressure 0.05-0.1 bar, which is more comparable to the industrial standards.

Figure 14.
Examples of more complex process configurations: 5-bed 5-step cycle with heavy reflux from the counter current depressurization on the left; 6-bed, 6-step cycle with feed, recovery step using the outlet stream of heavy reflux, heavy reflux, heavy reflux with light reflux product, counter current depressurization, light reflux and light product pressurization, on the right [290].
Alternatively, we can consider more complex multi-bed multi-step configurations. Below we review several studies that explore more complex process configurations in the context of post-combustion carbon capture. In particular, Reynolds et al. [291]., studied the capture of CO2 from a flue gas mixture containing 15% CO2, 10% H2O and rest N2 using potassium-hydrotalcite as the adsorbent. They had studied 9 different cycles with heavy and light reflux steps in 4-bed, 5-bed and 6-bed configurations.
Two examples of such advanced PSA/VSA processes are shown in Figure 14. A parametric study was then carried out and the best performing cycle was the 5-bed 5-step cycle with heavy reflux from the counter current depressurization (shown in Figure 14). The purity and recovery values were 98.7% with a productivity of 0.11 mol/m 3 .s.The next best cycle was cycle shown in Figure 14 on the right, which, although it showed a much lower recovery of 71%, had a high throughput of 1.11 mol/m 3 .s. It should be noted that this was a parametric study which did not show the optimum performance of these cycles, in other words, detailed process optimization was not carried out to identify conditions corresponding to maximum productivity while meeting purity and recovery targets.
Zhang and Webley [292] constructed a pyramidal hierarchy of cycles. The pyramid consisted of cycles ranging from a simple 2-step cycle to complex cycles which included heavy reflux, light reflux, pressure equalization etc. They carried out a parametric study on the effect of feed step duration, pressure equalization, rinse, evacuation and purge step on the purity and recovery values. Although, their model was a simpler one, nevertheless it provided some useful insights on the performance. They had studied the effect of feed time, light reflux, pressure equalization and heavy reflux steps on the purity, recovery and specific energy consumption. One of the main conclusions was that the addition of heavy reflux improved the purity from 85.7 % to 95.2% in their experiments.
A two-stage pressure swing adsorption process was studied with activated carbon by Shen et al. [293]. The first stage employed a Skarstrom cycle comprising of the following steps: pressurization with feed, adsorption, counter-current evacuation and light product purge. In the second stage a 5-step cycle comprising of pressurization with feed, adsorption, depressurizing pressure equalization, countercurrent evacuation and pressurizing pressure equalization steps, was used. The first stage used a feed under ambient pressure while the second stage required a compression up to 3.5 bar. With a vacuum pressure of 0.1 bar, the two stage process produced a 95% pure CO2 product with 74.4% CO2 recovery. Through a combination of experiments and simulations Wang et al. [294], studied CO2 capture from a dry flue gas containing 15-17% CO2. They used a two-stage process with the 1st stage containing 3 columns packed with zeolite 13X and the second stage containing 2 columns packed with activated carbon. In the first stage, the cycle chosen was an 8-step cycle comprising of pressurization with feed, adsorption, co-current evacuation, heavy reflux, depressurizing pressure equalization, countercurrent evacuation, light reflux and pressurizing pressure equalization. The CO2 product was collected from the counter-current evacuation and the reflux step. For the second stage, a six-step cycle comprising of pressurization with feed, adsorption, rinse, depressurizing pressure equalization, light reflux and pressurizing pressure equalization was used. The vacuum pressures of the first and second stage were 0.08 bar and 0.2 bar respectively. In the 1 st stage, CO2 purity of 70% was achieved. This stream containing 70% CO2 was then compressed and sent to the second stage and here the product contained over 95% CO2. The overall recovery was over 90%. From the experiments and the simulations, the values of the energy consumption were found to be 2.44 and 0.76 MJ/kg respectively. The large differences in the reported energy consumption values could have been a consequence of the low pressure used in the first stage which would have resulted in a lower vacuum pump efficiency as shown by Krishnamurthy et al. [80].
Haghpanah et al. [279] also performed detailed process optimization on 7-different cycles to identify the optimal configuration of post-combustion CO2 capture using zeolite 13X as adsorbent, ranging from 4 to 6 steps.
The genetic algorithm-based optimization was carried out in two steps: a) to maximise purity and recovery; b) minimize energy consumption and maximize productivity for cycles satisfying 90% purity and recovery constraints. The decision variables were the step durations, evacuation pressures and the feed flow rate. The adsorption step pressure was kept constant at 1 bar and the feed was a 15% CO2 and 85% N2 mixture at 298K. The optimization results from the 1 st step showed that 4-cycles namely 4-step cycle with LPP, 5 step cycle with light reflux (LR) and LPP and the two 6-step cycles satisfied the 90% purity-recovery targets. The next step was to minimize energy and maximize productivity and, in this case, the 4-step cycle with LPP was the best performing cycle in terms of the energy consumption. The minimum energy consumption was 0.47 MJ/kg CO2 for a productivity of 0.57 mol/m 3 .s and with an evacuation pressure of 0.03 bar. The 6-step cycle with LR and HR was the best in terms of the purity and recovery and in this cycle over 98% purity and recovery were achieved. However, this cycle had much higher energy consumption. The 4-step cycle with LPP was also able to achieve the 95% purity and 90% recovery targets and in this case the energy and the productivity values were 0.554 MJ/kg CO2 and 0.44 mol/m 3 .s respectively. The 4-step cycle with LPP was shown to meet the 95% purity and 90% recovery targets through a pilot plant study by Krishnamurthy et al [80].
In a later study, Haghpanah et al. [207] studied CO2 capture using a carbon molecular sieve by the optimization of 1-stage and 2-stage VSA processes. The two-stage process basically is the 4-step cycle with LPP carried out twice with the product from the counter-current evacuation step serving as the feed for the adsorption step in the second stage. In the one stage process, 5-step cycles with heavy reflux and with feed and light product pressurization were studied. It was seen that the 2-stage VSA process was the best in terms of energy and productivity. However, the productivity was about 50% smaller than that of zeolite 13X mentioned above as carbon molecular sieve is a kinetically selective adsorbent.
It should be noted that in both the studies by Haghpanah et al. [207,279] the evacuation pressures were from 0.01 to 0.05 bar and this, as we have already discussed, may not be industrially achievable. In a recent study by Khurana and Farooq [284], it was shown that with a 6-step cycle it was possible to achieve evacuation pressures of 0.1 bar and above. The 6-step cycle is essentially the 5-bed 5-step cycle with heavy reflux from light reflux product of Reynolds et al., and with the addition of a cocurrent evacuation step [291]. The authors have compared the performance of this cycle with that of the 4-step cycle with LPP and used two adsorbents namely zeolite 13X and UTSA-16. Through detailed process optimization, it was seen that the 6-step cycle was able to achieve similar productivity values as the four-step cycle but at a much higher evacuation pressure of 0.1 bar.
Another cycle that has shown promise for producing both CO2 and N2 in high purities is the dual reflux pressure swing adsorption cycle. The concept of the Dual reflux PSA (DRPSA) process first appeared in the 90's in the work by Diagne et al. [295,296], who had obtained 95% CO2 purity and recovery from a stream containing 20% CO2. The DRPSA contains two columns one operating at high pressure and the other at low pressure at a given instance. Feed can be introduced from the bottom or from an intermediate position and both at low and high pressures. The enriched gas from low pressure feed is then compressed and sent to the column at high pressure to perform the heavy reflux. The light product from this heavy reflux step can be used to recover the heavy component simultaneously during the feed step. After this step, the column roles reverse, and the same sequence of steps are carried out. Over the years variations of the dual reflux PSA cycle have been studied by a few authors [297][298][299][300]. One among them is that of Li et al. [299], who had studied the CO2 capture from a binary mixture containing 15% CO2 and 85% N2 at 2 bar and 20°C feed using silica gel adsorbent. They were able to achieve over 99% purity and recovery for CO2 and N2. In a follow up work Shen et al. [300], used the same cycle and the adsorbent and achieved over 95% purity and recovery with an energy consumption of 2.5 MJ/Kg.
In Table 8 below we complement the summary of Abanades et al. [267], with the summary of the recent studies of various process configurations for post-combustion capture. We note that although contains two TSA cycles, in reality the productivity of these cycles will be low compared to the PSA/VSA processes: this would lead to columns simply too large (or to a large number of them) to be competitive in the post-combustion capture from coal-fired plants. They may however find application in the carbon capture from natural gas fired power plants, where the concentration of CO2 in the flue gas is much lower. Naturally, a question emerges what process and material we can already identify as the most promising for post-combustion carbon capture. As it turns out, this is a challenging question simply because previous studies have used different process and cycle configurations for assessment of materials performance so that consistent comparison of the results produced by these studies is not possible.
For this reason, we constrain ourselves specifically to the 4-step VSA-LPP process configuration with the following conditions of the feed: composition: 15% CO2 and 85% N2, temperature = 298.15 K, pressure = 1 atmosphere. For this, we also impose the purity and recovery requirements to be equal to 95% and 90% for CO2 respectively. For a long time, zeolite 13X remained the key benchmark material for this process with the lowest energy consumption reported at 0.53 MJ/kg of CO2 captured [78]. Using the same process, Khurana and Farooq obtained 0.45 MJ/kg CO2 captured for UTSA-16 MOF materials [324]. Using the same model as a benchmark, Balashankar [82] identified several hypothetical zeolites and ZIFs, performing better than 13X (Figure 15.). A slightly different optimization target was used by Burns et al. [84], where a similar multiscale screening strategy was applied to screen more than 5000 MOFs from the Core-MOF database (using the heuristic filters 1584 MOFs from the database were passed on to the further investigation using process modelling tools). In this case the target was to identify all materials that simultaneously achieve a parasitic energy less than 0.903 MJ/kg CO2 (this benchmark corresponds to the solvent based capture energy penalty) and a productivity greater than that of zeolite 13X (4. NB. All energy values are electric, th: thermal, mix: electric+thermal. CO2/m 3 ). Indeed, as seen from Table 9, several materials are able to achieve this target including UTSA-16 which is also previously identified as a promising candidate [50]. Below, in Table 9 we provide the data from " Table 1" of the original publication in units consistent with Table 8 and throughout the article.

Emerging Numerical Techniques for Process Optimization
The process simulations we have covered so far are computationally expensive: a single process simulation for a given set of design variables takes minutes to complete. Process optimization to obtain a Pareto front as described in Figure 14 requires thousands and tens of thousands of simulations, leading to an overall cost of the process optimization exercise of 10 2 -10 3 CPU hours for a single material. Clearly, routine screening of tens, hundreds or thousands of materials at the process level is prohibitive.
This promoted the development of several strategies to reduce this cost. These strategies can be split into three main categories: 1. Reduce the pool of candidate materials by low cost, preliminary screening strategies 2. Reduce the computational complexity of the individual process simulations in the optimisation process: a. Accelerate the convergence to CSS b. Use a simpler model from the model hierarchy c. Replace the high-fidelity model with a surrogate model trained on the high-fidelity model 3. Reduce the computational effort of the optimisation process The three approaches can be combined, but all have disadvantages and limitations, which can compromise the screening process so that the optimal material and cycle configuration is missed. Here we review studies with a focus on accelerating process optimization using strategies outlined above.
Strategies in the first category can use any of the previously published performance metrics to eliminate candidates from the candidate pool so that the expensive computational effort can be spent on the most promising candidates. As described in Section 4, simple performance metrics are not able to correctly and accurately rank materials for the complex and highly dynamic adsorption processes where the performance is given by a balance between the competing objectives of energy penalty and productivity as well as the competing constraints of purity and recovery. Thus, it is crucial to have very conservative exclusion criteria so that potentially promising candidates are not removed from the candidate pool. On the other hand, a number of these metrics can be computed very quickly so that the least promising candidates can be removed for a low computational cost.
Burns et al. [84]. performed a detailed multi-objective process optimisation and ranking for a large range of materials for post-combustion capture. Afterwards, they trained Machine Learning (ML) classifiers to predict the objectives, i.e. purity, recovery, parasitic energy and productivity, based on 29 sorbent metrics such as working capacity, selectivity, and isotherm parameters. They showed that the N2 adsorption behaviour is crucial for the correct classification of materials that meet the 95/90 purity/recovery constraints and achieved a prediction accuracy of 91% for this. However, the prediction of parasitic energy and productivity for materials that achieve the 95/90 purity/recovery constraints achieved only very low prediction accuracy. They concluded that full process simulations are required for accurate ranking of parasitic energy and productivity. An interesting approach is followed by Khurana and Farooq [50] who trained a classification neural network based on five equilibrium isotherm characteristics which cover the parameter space of the dual-site Langmuir isotherm. Their model can predict with 94% accuracy whether a material can meet the 95/90 purity/recovery constraints for post-combustion capture with the LPP-VSA cycle. In addition to the good accuracy, the false negatives, i.e. materials that can achieve the target but were wrongly classified, showed high energy penalties and low productivity. For the materials that met the 95/90 purity/recovery constraints, they developed a meta-model to predict the energy penalty and productivity and achieved R 2 values of around 0.9 for minimum energy penalty and maximum productivity.
The second category is split into three methods to reduce the computational complexity of the individual process simulations. First, instead of simulating cycle after cycle to reach CSS, so-called successive substitution, several studies have explored methods to accelerate the convergence to CSS. For example, Smith and Westerberg [328] and Ding and LeVan [329] used Newton and quasi-Newton steps to reduce the cyclic deviation. This method requires the calculation of the Jacobian and can achieve about an order of magnitude faster convergence. Alternatively, derivative-free extrapolation methods such as the epsilon extrapolation used by Friedrich et al. [330] can reduce the required number of cycles to CSS by a factor of 3. Pai et al. [331] use artificial neural networks to predict the bed profiles at CSS and use this to initialise the high-fidelity simulations. In their tests it reduces the average number of cycles which need to be simulated to reach CSS by a factor of 6.
Second, simpler but still physics-based models are used instead of the high-fidelity models. These simplified models should be fast to calculate while still capturing the main physics of the separation process. Balashankar et al., [332] use a batch adsorber analogue model as a simplification for the full VSA model with spatial discretisation. The simplified model assumes that the system is isothermal, well-mixed and has no mass transfer resistance, but still captures part of the physics of the separation and can be solved in seconds. The authors compared the output from the simplified model with the detailed process optimisations and developed a classifier that achieved a Matthew correlation coefficient of 0.76 in the classification of materials that meet the 95/90 purity/recovery constraints. In addition, they calculated a linear regression for the energy penalty, which estimated the energy penalty with reasonably good accuracy, i.e. within 15% for 83% of the materials. However, Biegler et al. [333] evaluated the use of simplified models for process optimisation and concluded that it can lead to convergence failure and even to false optima.
Third, surrogate models are built based on the output of the high-fidelity models. These models are faster to evaluate and are usually embedded into optimisation methods. Agarwal et al. [14] used proper orthogonal decomposition (POD) to replace the detailed spatial discretisation with a reduced order model (ROM), leading to a system of Differential Algebraic Equations (DAE) of a significantly lower order. The ROM was trained on a number of bed profiles for different cycle conditions simulated to CSS. Because only the largest singular values are used for the ROM, the size of the discretised model is reduced by an order of magnitude. The ROM is accurate close to the training cycle conditions but loses accuracy further away from these points. This means that the ROM needs to be retrained if the optimisation moves away from the original training points.
In recent years, the focus has moved to directly using the optimisation objectives and constraints to build ROM instead of using ROM of the bed profiles. This approach replaces the process simulation (reduced order or high-fidelity) with fast-to-calculate surrogate models (also called ROM, meta models or emulators), which directly calculate the optimisation objectives based on the optimisation variables. These models are built from the input-output relations generated with the high-fidelity models and can be used with any black-box optimisation algorithm. This enables the interfacing with state-of-the-art multi-objective optimisation methods to handle the trade-off between competing objectives and constraints.
The process of a surrogate optimisation is shown in Figure 16. The process starts with an initial Design of Experiments (DoE) which should cover the entire design space. The high-fidelity model is used to simulate the responses for these initial designs. Then the optimisation loop starts by building a surrogate model based on these input-output relations. The optimisation method operates on this, fast-to-calculate, surrogate model to find promising design points. The choice of the next design points is a balance between exploring the design space and exploiting the best-predicted design or designs. The new design point is evaluated with the high-fidelity model, added to the input-output relations and a new iteration of the optimisation loop starts, i.e. we build a new surrogate model. The optimisation loop is stopped once a stopping criterion, which is often a computational budget, is fulfilled. Beck et al. [15,265] used Kriging regression based surrogate models with the NSGA-II optimiser to simultaneously optimise the CO2 purity and recovery for post-combustion capture. The Kriging regression models the input-output relation as a Gaussian process and gives the best linear unbiased prediction. In addition, it also provides confidence bands for the prediction, which can be used to explore the design space. They achieved a five-times reduction in computational effort and also investigated the specific energy penalty [265].
The rapid development of machine learning methods and, in particular, artificial neural networks (ANN) is mirrored in the application of machine learning to adsorption process optimisations. Sant Anna et al. [334] developed a three layer neural network (input layer, one hidden layer, and output layer) surrogate model for the separation of nitrogen and methane. They trained the neural network on around 500 training samples and performed a multi-objective optimisation of N2 purity and recovery on the trained network without further updating the surrogate model. Comparing the optimal values with the high-fidelity simulations showed that the maximum relative difference was 1.4% for N2 purity and 4% for N2 recovery.
Instead of directly approximating the optimisation objectives and constraints, Leperi et al. [335] used ANN based surrogate models to approximate each basic step, e.g. counter-current pressurisation and co-current feed, of the PSA cycles. This approach enabled them to build arbitrary PSA cycles and to include cycle synthesis in the optimisation procedure without the need to retrain the ANN for each process configuration. They built 12 surrogate models for each step: one ANN for the state variables at 10 locations along the column and one for each end of the column to predict the inflow/outflows during the step. The ANNs are trained with high-fidelity simulations for 300 Latin hypercube samples and used to predict the column profiles as well as purity and recovery for three process configurations and two adsorbents for post-combustion capture. The predictions were used in an optimisation loop to find the purity/recovery Pareto front. The solutions on the Pareto front were used to test the accuracy of the ANN prediction and to retrain the ANN in case the prediction is too far from the highfidelity simulation. After retraining, the relative errors for both purity and recovery were below 1.5% for all cases.
Subraveti et al. [336] used a surrogate model based on ANN in the multi-objective optimisation of purity and recovery of pre-combustion CO2 separation and achieved a 10 fold reduction in computational effort. For the first five generations of the multi-objective NSGA-II algorithm they used the high-fidelity model. This generated training data for the ANN, which would already be biased towards the optimal region of the design space and should improve the prediction accuracy in the optimal region. The ANN with one hidden layer with 10 neurons was trained on the generated inputoutput data. The remaining 45 generations of the optimisation were performed on the ANN. The Pareto front was close to the one generated with the high-fidelity model, but had a relative error around 1% in both objectives. In a subsequent paper, the group compared a range of machine learning methods and showed that Gaussian process regression achieves an R 2 value above 0.98 for purity, recovery, energy penalty and productivity with a training set of 400 randomly sampled high-fidelity simulations [331]. Their optimisation on this surrogate model (without further refinement) was within 3% of the high-fidelity simulation for purity and recovery as well as for energy penalty and productivity. However, the latter was for a reduced 95/80 purity/recovery constraints.
Pai et al. [86] developed a material-agnostic surrogate model called MAPLE that fully emulates operation of the 4 steps VSA-LPP cycle at the cyclic steady state. The framework is based on a dense feedforward neural network trained with a Bayesian regularization technique. The framework accepts the adsorbent properties, the Langmuir adsorption isotherm parameters, and operating conditions as input. It predicts key performance indicators of the process including CO2 purity and recovery in addition to productivity and overall energy consumption of the process as output. The model was trained by a set of data generated using detailed process modelling. In order to reduce computational time of the multi-objective optimization, MAPLE was used to calculate the CSS performance indicators and feed them back to the optimizer. The fully trained model, predicts each performance indicators with less than 2% error compared to the detailed process modelling. The computational time required for simulation and optimization of the process was also reduced from 1500 core-hours per adsorbent to ≤1 core-min for each adsorbent which shows a significant improvement for screening of large databases of porous materials [86].
Strategies in the third category include a range of methods to reduce the computational effort of the optimisation method itself, i.e. reduce the number of required iterations to reach an optimal value or Pareto front. The first strategy should be the reduction of the search space. This includes the removal of parameters, which have no or only a small impact on the performance and the reduction of design space, i.e. reduce the evacuation pressure range. For example, Balashankar et al. [332] removed the blowdown and evacuation times from the list of optimisation variables. This was acceptable in their optimisation because these variables have very limited impact on the purity and recovery. However, they have a large effect on productivity and energy penalty.
Yancy-Caballero et al. [85] performed a hierarchical, multi-objective optimisation with NSGA-II. They first optimised purity and recovery to screen for materials that achieve the 95/90 purity/recovery constraints and then optimised the promising materials for energy penalty and productivity. The energy penalty and productivity optimisation was seeded with the results from the initial optimisation and was performed in two steps: the first step used a low spatial resolution, which reduces the computational complexity, and the second step used a high spatial resolution and was pre-seeded with the low resolution results.
Finally, Ding et al. [337] and Jiang et al. [338] presented a strategy which combines the reduction of the computational complexity of individual simulations with a reduction of the computational complexity of the optimisation. They included the CSS condition as a constraint in the constrained single-objective optimisation problem so that both the objective and approach to CSS were optimised simultaneously. This approach, called the simultaneous tailored approach by Biegler and co-workers [7], removes the expensive calculation of CSS for each iteration and has reduced the computational time by a factor of 10 for single-objective optimisations of air separation VSA cycles [338].

Available Tools and Software for Process Modelling and Optimization
The objective of this section is to introduce the reader to several software packages and libraries that are available for PSA/VSA process simulations. Broadly these can be divided into three categories: codes developed by the academic groups, codes developed within various companies and commercial software packages, with build-in adsorption process simulators. From this classification we can also immediate identify the most significant challenge in a consistent description of these tools: the academic and industrial codes are in-house developments, used by the academic or industrial groups and we do not direct have access to the organization, functionality, implemented models or capabilities of these codes to make the comparison consistent. In the case of the commercial codes, although these are typically accompanied with a sufficient documentation and case studies, we do not have access to the organization of the codes and we do not believe a systematic comparison of the commercial packages in terms of computational efficiency, performance and scope has been carried out. Hence, from the onset we admit that this section is likely to be incomplete. We will first review the software packages, and then we will briefly touch on the platforms (programming languages, libraries, software within which adsorption process simulators can be developed.

Commercial software:
gPROMS: The process builder developed by Process Systems Enterprise (PSE) has an adsorption process library which has been used for simulation of pressure and temperature swing adsorption processes [313,327]. In the adsorption process library, it is possible to use the dispersed plug flow or the plug flow model. The adsorption isotherms of Langmuir, dual-site Langmuir and virial isotherms can be used in gPROMS. Here, the flow sheet can be built by joining individual units such as valves, headers mass flow controllers, sources and sinks and adsorption columns. The adsorption process model is a system of partial differential equations that are discretized in the spatial domain using either finite difference (backward, forward and central), finite element, finite volume or orthogonal collocation with finite element. The flow controllers supply constant amount of gas, while the sources and sinks are used to specify initial and final operating conditions. gPROMS also has the facility to account for column headers to distribute flow and these are modelled as continuous stirred tank reactors. Building a flow sheet enables one to schedule various steps operating in multiple columns. In principle, within gPROMS it is also possible to perform optimization and scheduling of VSA process using in house libraries [313,327]. PSA processes have been optimized in gPROMS for the maximization of CO2 product purity and recovery, with the number of beds, process configuration, feed pressure, particle diameter, length to diameter ratio and feed flowrate as the decision variables by Nikolaidis et al. [339]. The same approach was used in an earlier study by Nikolic et al. [340] for H2 recovery from steam methane reformer off gas. However, it should be noted that these studies have not reported any Pareto fronts.
Aspen Adsorption: It is a flow sheet simulator that can design, simulate and optimize adsorption processes [341]. Few studies exist in literature that have simulated PSA and dual-reflux PSA processes for CO2 capture using this program [297,299,300,322]. In Aspen Adsorption, it is possible to simulate multi-bed PSA processes with isothermal or non-isothermal model and use non-ideal gas equations of states. In most publications with Aspen Adsorption, a Langmuir model has been used. Moreover, it is also possible to use finite difference or finite volume numerical schemes to solve the model equations.
To the best of our knowledge no cycle optimization studies have been published using Aspen Adsorption software, although it is possible to couple Aspen products and MATLAB [342].
ProSim DAC: ProSim DAC is a dynamic simulation software from ProSim 3 . It is capable of simulating adsorption and desorption steps using TSA, PSA and VTSA processes. Data for a wide variety of adsorbents is available (e.g., activated carbon and zeolites) and is accompanied by many different models for equilibrium data (adsorption isotherms) and mass transfer models. As this is a relatively new addition to the ProSim family of process simulation tools, we are not aware of any academic article on carbon capture simulation and optimization that have used ProSim DAC.

Academic codes:
Several research groups (Yang, Ritter, Ebner, Mazzotti, Webley, Rodrigues) have been developing codes for adsorption process modelling starting from the 1980ties. The key challenge in the discussion of the codes used by the academic group is similar to the issues associated with the industrial software: the codes are usually not open source, and full details of the algorithms employed, implementation and capabilities are not available. Hence, it is impossible to comprehensively comment on these codes, not to mention to compare them on a consistent basis. Here we mention some of the codes we are aware of: MINSA (Monash Integrated Numerical Simulator for Adsorption): MINSA is a generalized cycle simulator which has been developed by Webley and co-workers for PSA simulations using the DVODE integration scheme of Brown, Bryne and Hindmarsh [343] written in FORTRAN [307,[344][345][346]. This simulation package solves mass and energy balance equations that have been discretized by the finite volume method [307,346]. The software has been used extensively for various adsorption processes and verified against experimental data over the past two decades [292,345,347,348] CySim (Cycle Simulator): CySim is a modular computer program for simulation of adsorption processes which has been developed by Brandani and co-workers [330,349] in University of Edinburgh. CySim can be used to simulate breakthrough curves, ZLC experiments [349], dual piston PSA and other PSA processes. The user defined structure is translated into a system of differential algebraic equations, which are solved with the SUNDIALS library. This can be interfaced with either MATLAB or Python's genetic algorithm packages such as gamultiobj [350,351], inspyred [352] and Platypus [286] to perform process optimization. Recently Farmahini et al. have used CySim to simulate and optimize the 4-step VSA process with LPP for post-combustion carbon capture [51,83]. CySim is regularly updated with new models and applications for example for monolithic adsorbents to include inlet and flow maldistributions [283,353].

Code by Yancy-Caballero et al.:
Recently, Yancy-Caballero and co-workers developed a MatLab code to simulate PSA/VSA processes [85]. In particular, the code uses the finite volume method with the weighted essentially nonoscillatory (WENO) scheme to discretize the PSA model; and the ode15s solver within MATLAB to solve the resulting ODEs. NSGA-II algorithms within MATLAB is employed for process optimization. In the most recent study, this code has been applied for performance ranking of several MOF materials in post-combustion carbon capture processes, using Skarstrom, a fractionated vacuum swing adsorption (FVSA), and a five-step cycles. A notable feature of the code is that it is open source and is publically available from the github depository of the academic group which has developed the code. Table 10 provides a list of process simulation software that have been discussed in this section.

Existing Challenges and Open Questions in Multiscale Screening of Porous Materials
In this section, we outline what we believe are key challenges in the implementation of the multiscale workflows for adsorption applications. Our awareness of these challenges evolved over time, as we being a collaborative group of molecular simulators and process simulators ourselves, navigated this emerging research field

Data Availability and Consistent Implementation of Multiscale Screening Workflows
Imagine, two articles have been published by two different academic groups ranking two sets of porous materials for a particular separation using the multiscale workflows discussed in this review. Are these rankings compatible and consistent with each other? In other words, can the results of the two studies be combined in one ranking, and therefore the best material identified out of the combined group of materials in both articles? In general, the answer is no. From the studies we have reviewed in this article, it is clear that different groups employ different hierarchies of models and assumptions, use different sources and values of parameters, and apply different conditions of the process. As we have discussed in Section 6.3.5, consistent comparison is possible only on a limited subset of available studies. Here in this section, we would like explore the issue of data availability and consistency of simulations in more detail.
As has been discussed throughout the article, not all of the data required to set up multiscale process simulations is available from molecular simulations. This is in fact a limiting factor for multiscale materials screenings to be performed fully in silico. Some of the data required for process modelling such as, for example, diffusion models for the bulk gas phase, can be constructed using appropriate, well-established theories and models. For other data however, some assumptions have to be made.
Recent studies have probed the influence of some of these assumptions on the actual process performance predictions as outlined below: Heat capacity of the adsorbent: Traditionally, it has been assumed that heat effects play a secondary role in the adsorption column and will not significantly influence the performance of the process or ranking of the material. Adsorbent heat capacity is also a property scarcely measured or available for the new classes of porous materials, such as MOFs. As a result, a pragmatic approach has been to assume the heat capacity of all materials to be equal to the heat capacity of a reference material, such as Zeolite 13X [50,77]. Recent preliminary sensitivity study of the influence of this parameter, painted a somewhat different picture: performance of HKUST with the value of the heat capacity equal to that of zeolite 13X was considerably different from the performance of the same material, using an experimental value of this property [83]. So it seems for a diverse group of porous materials (MOFs, zeolites), it would be prudent to procure more heat capacity data, explore these heat effects in more detail and use more accurate data in multiscale screening studies. How can this be accomplished in purely in silico workflows? In principle, we can use computational phonon analysis methods [254,255,357,358] to generate heat capacity of individual materials. This, however, will add yet another level of complexity to the already long chain of simulations that must be performed to compute the prerequisite properties for multiscale simulation of PSA/VSA systems. Analysis of phonon properties of porous solids could require time-consuming quantum mechanical calculations, considering that currently available classical force fields may not necessarily be able to capture the full range of harmonic/anharmonic properties pertinent to thermal properties of new classes of materials [358,359]. Therefore, the remaining challenge here is to develop affordable computational techniques (such as group contribution methods [360]) by which thermal properties of porous materials can be calculated without compromising the accuracy of the predictions.
Pellet size and pellet porosity: In the traditional process modelling literature, the values of pellet size and pellet porosity are typically obtained from experiments for a specific sample of a material under consideration. The values of these characteristics cannot be provided by molecular simulations or from some approximate theories. Again, the pragmatic approach, adopted in the previous studies, has been to use the values available for a reference material, such as Zeolite 13X. However, in the context of the ranking of porous materials, a question can be asked whether optimal performance of a material can be achieved only at material-specific values of pellet size and pellet porosity (within the range of feasible experimental values)? In this case, shall the ranking be performed under the constraint of specific values of pellet size and porosity, or shall these properties also become some optimization parameters? Farmahini et al. [83] have recently explored these questions in a recent study [83], and observed that depending on the protocol (pellet size and porosity are constrained to the reference values of Zeolite 13X versus being free optimization parameters) the order of top performing materials does change [51,83]. Therefore, consistent comparison of the materials' performance at the process level must take into account opportunities for the geometry optimization of the column packing.
Process level models for adsorption isotherms: Adsorption models for process simulations are trained on the available experimental or simulation data and they should be able to give consistent and accurate representation of multicomponent adsorption equilibria. This is however not always the case and two different models trained on the same data may give different predictions of the binary (or multicomponent) adsorption equilibria (depending on the training protocols adopted) and hence, process level predictions [51,361]. For example, Farmahini et al. [51] have recently demonstrated that the use of different recipes for fitting adsorption data to the DSL model can affect position of the energy-productivity Pareto fronts obtained from the process optimization. The authors have therefore proposed and validated a rigorous numerical protocol for consistent fitting of the widely used DSL model, in which the choice of temperature range, fitting constraints, and calculation of Henry's constant are standardized [51]. Similarly, several other studies have attempted to establish consistent routines for fitting equilibrium adsorption data [361][362][363], nevertheless none of the proposed procedures have been universally adopted by other groups, as a result of which consistency of various materials rankings that have been produced so far remains uncertain. As has been also discussed by Farmahini et al. [51], the ultimate test of the analytical models used in the process level simulations is their ability to predict binary and multicomponent adsorption data. This can be easily done using molecular simulations, as simulations of multicomponent systems come with relatively small additional effort compared to the simulation of single component systems. Hence, we propose that this validation step is used to probe the accuracy of the analytical models for adsorption isotherms before performing any process simulation.

Sensitivity analysis and propagation of errors:
From the studies reviewed so far, it is clear that overall process performance and ranking of materials depends upon calculation of a large group of parameters and model assumptions at both molecular and process level of descriptions (see Table 6). Despite some studies on the sensitivity of process performance to various input parameters and data [51,83,85,239,364], it is yet to be established what level of accuracy is required for the full spectrum of parameters and models to guarantee consistent and comparable ranking of porous materials between different studies. One crucial element of such a study would be the investigation of errors propagation from molecular level all the way through to process modelling and optimization. For example, we can understand how the errors arising from the use of inaccurate atomic force fields in GCMC simulations for prediction of adsorption isotherms are combined with the errors resulting from the use of numerical models for fitting adsorption data, and what impact they will have on the overall performance of the process.
Consistency between various simulation packages: One important aspect in developing consistent multiscale workflows for materials ranking is having access to widely-used and validated open-source software and case studies. As mentioned in Sections 6.1.2, 6.2.2 and 6.2.4, there are a number opensource molecular simulation packages for which several benchmark and case studies have been produced [145,159,208]. However, this has not been done for process simulation packages mainly because the majority of these software are not available as open-source. In fact, among all the software introduced in Section 6.3.7, only one code has been released as open source [85]. As shown by the molecular simulation community, having access to open-source software and clear case studies will facilitate expedited development of new generations of software and tools for material screening studies. We believe this can be also true for the process simulation community and hope that the current review has been successful in demonstrating the importance of any efforts which can address the current gap.

Improving Efficiency of Process Optimization for Comprehensive Screening of Materials Space
Multiscale simulation of PSA/VSA processes for screening of large databases of porous materials requires extensive computational resources. In the screening workflow, process optimization is usually considered as a bottleneck where significant computational efforts are incurred [336]. Attempts have been made to improve computational efficiency of process optimization, through reducing dimensionality of the variable space in process optimization, by development of novel machine-learning methods (Section 6.3.6), or by hierarchical approaches, where the simplified process models are used first in the pre-screening studies, while the full process optimization focuses on a smaller subset of cases. Together, these methods pave the way not only for faster screening of large databases of porous materials but also for identifying the most efficient process configurations for a particular separation process. This can lead to better understanding of the material-processperformance relationships. Nevertheless, the remaining challenge is yet to tackle the magnitude of the material-process phase space. Currently, these methods have only been tested for screening of small sets of porous materials (<2,000) [84] which is infinitesimal compared to the huge number of materials that has been discovered so far, as mentioned in Section 6.1.1. Also, experimental evidence for validation of numerical techniques that are used for expedited optimization of PSA/VSA processes are still scarce [331] and it is for the future studies to address this important limitation.

Prediction of Equilibrium Adsorption Data for Nitrogen
From the perspective of a process modeller, the most immediate advantage of having access to accurate molecular simulation tools is gained by having the ability to predict multicomponent adsorption equilibria. This requires force fields for the mixture constituents that are validated against pure component adsorption data of good quality. Accurate adsorption equilibrium models are the basis for equilibrium driven separations, such as those considered here for post combustion capture.
As mentioned in earlier sections of this review, several studies have recently proposed that accurate prediction of nitrogen adsorption data plays an important role in correct prediction of material performance at the process level [6,50,51,55,79,84,239]. For example, these studies indicate that it is the nitrogen adsorption behaviour that most significantly determines whether the process can produce CO2 with 95% purity and 90% recovery or not [6,55,84]. Even for those materials that meet 95%-90% criterion for purity and recovery of CO2, inaccurate prediction of N2 adsorption data was shown to affect energy and productivity of the process [51]. These observations in turn highlight the importance of having access to accurate atomic force fields for more accurate prediction of nitrogen adsorption data using molecular simulations, the topic whose importance has been overlooked so far. Force field development for nitrogen is, however, a challenging task for two related reasons. Initially, one would consider QM methods to develop a detailed picture on the potential energy surface for nitrogen in various porous materials including MOFs, in a similar fashion that has been done for several CO2-MOF systems. What is important to remember is that nitrogen is a weakly adsorbing component (heat of adsorption 10-20 KJ/mol, but likely to be closer to 10-12 kJ/mol). Relative error in QM estimates of energy of binding is likely to have a much stronger impact on nitrogen adsorption than on stronger adsorbing carbon dioxide. For a similar reason, the uncertainty in the experimental adsorption isotherms for nitrogen is also greater as the amount adsorbed tends to be small under conditions of interest. As a result, and as discussed in Section 6.2.5, the number of QM-derived force fields which are especially developed for accurate prediction of nitrogen adsorption in novel classes of porous materials is extremely limited [240][241][242][243]. At the same time, there are numerous examples in the literature suggesting inadequacy of generic force fields for accurate prediction of experimental nitrogen adsorption in many materials including zeolites, MOFs and ZIFs [51,83,240,[244][245][246]365]. Given the importance of this topic for accurate prediction of materials performance at the process level, some strategies for the development of reliable force fields for nitrogen adsorption need to be proposed.

Validation of Multiscale Screening Workflows
Despite recent advances in development of more sophisticated multiscale screening workflows, validation of the materials rankings produced by these frameworks is still an outstanding issue. As illustrated throughout this review, multiscale screening workflows have a modular structure in which various computational modules are put together to perform different types of simulations. The simplest workflow contains three modules in order to perform (1) GCMC simulation, (2) process modelling, and (3) process optimization. This can be further extended, if one decides to include quantum mechanical or MD simulations in the workflow. Normally, results of each module can be validated separately. For example, adsorption isotherms generated using GCMC simulation are routinely compared against equilibrium adsorption data obtained from experiments to ensure the accuracy of the predictions. At the process level, validation tests are conventionally carried out by reproducing column breakthrough curves, temperature, pressure and concentration profiles from experiments [80,81,330,366,367]. Efforts for validation of genetic algorithms which are used for multi-objective optimization of PSA/VSA processes are more recent and less wide-spread. In one recent example, the ability of multi-objective optimization techniques to guide the design of PSA/VSA processes have been shown by Perez et al. [81]. In this study, purity and recovery of CO2 predicted through numerical optimization of a basic 4-step VSA cycle and a 4-step VSA cycle with LPP were replicated by experiment for post-combustion carbon capture using zeolite 13X as adsorbent [81]. Unfortunately, it was not possible in this study to carry out any measurement to verify total energy consumption of the process against experimental data [81]. A pilot plant study conducted by Krishnamurthy et al. [80] for CO2 capture using the same 4-step processes with zeolite 13X report significant discrepancies between theoretical and experimental values of energy consumptions, while their analyses show overall quantitative agreement for purity and recovery, and somewhat modest agreement for productivity data [80]. Other pilot-or lab-scale studies also report reasonable agreement between experimental and theoretical values of purity and recovery at the given feed concentrations for separations of CO2 using VSA/PSA processes [317,318], while the discrepancy reported for the total energy consumption is still considerable [318]. From the above studies, it is clear that full validation of process optimization module against reliable experimental data is still an open issue particularly with regard to correct prediction of total energy consumption of the process. This becomes especially crucial, if we remember that energy-productivity Pareto fronts obtained from multi-objective optimization of the process play a central role in performance-based ranking of porous materials. With the surge in development of machine-learning approaches for modelling and optimization of separation processes, one would also need to consider additional tests for validation of these novel techniques. Some recent studies have reported promising cases where machinelearning based surrogate models developed for the optimization of PSA/VSA processes have been validated against Pareto fronts and cyclic steady state (CSS) column profiles that are mainly obtained from detailed process simulations [86,331,335] but some also from lab-scale experiments [331].
Despite all these efforts for validation of different computational modules of multiscale screening workflows, there is no single material ranking study in which the order of top performing materials has been confirmed experimentally. In fact, unless this final level of validation is achieved, it is unlikely that the top-performing materials proposed by various computational screening studies are found their way into any industrial application.

Development of More Advanced Workflows
In addition to what has been discussed in this section, development of more advanced multiscale workflows for PSA/VSA/TSA processes can be envisioned where behaviour of more complex materials is simulated. An important example of such cases is the prediction of separation performance of many novel porous materials [247,248,368] such as materials with gaiting effetcs and phase-change adsorbents that exhibit step-shaped adsorption isotherms [270,[369][370][371]. Atomistic structures of these materials undergo considerable structural changes in response to external stimuli such as heat, pressure, humidity and adsorption of guest molecules [157,369]. Simulation of adsorption process in this class of materials must capture the interplay between presence of adsorbate molecules and structural deformation of the framework using computational methods that go beyond conventional GCMC (e.g. the osmotic Monte Carlo method or hybrid MC/MD methods) [144,372]. In addition to simulation of structural flexibility of these materials that must be handled at the molecular level, it is also crucial to develop more sophisticated analytical adsorption models that can capture stepwise shape of the isotherms in these materials as required for process simulation [270,373,374]. These two issues alone pose a significant challenge to the development of the future generation of multiscale simulation workflows for screening of flexible materials.

Future Outlook and Final Remarks
In this article we reviewed the recent progress in the application of multiscale workflows for material screening in post-combustion carbon capture. To make it useful for a wide audience of material scientists, chemical engineers and computational modellers, we introduced the basic principles involved in each element of the workflow and provided references to the available computational tools. We outlined what data is obtained at each level and what data is required at the next level, including the issue of availability and completeness of the data. The article also summarized all the recent studies in the field, and as such can serve as starting point for further developments. Reflecting on the state-of-the-art in the field, our take home messages are as follows: Beyond post-combustion carbon capture: In this article we focused on the post-combustion carbon capture as it is a very challenging, societally relevant and most investigated process. However, we believe the multiscale screening approaches reviewed here will become a new way to design and appraise material options in other separation applications. Specifically, when considering the need to decarbonize the chemical industry by 2050, it is clear that many new carbon capture materials and processes will need to be developed. The key aspects to consider are application of a wider range process configurations, and process conditions (primarily different levels of carbon dioxide concentration), and the relatively small scale of the processes, compared to post-combustion capture for energy generation (meaning, smaller amounts of the material are required). For these processes, it is likely that faster cycles will be used to reduce the footprint of the units, especially in retrofit applications. This will in turn lead to larger effects of mass transfer and heat transfer limitationsprecisely the challenges that need to be explored within the multiscale framework.
Air separation is a very useful case to consider for benchmarking multiscale modelling approaches. Production of oxygen is an equilibrium driven separation where the light component is produced. Therefore, simpler process configurations will work well in this case and advances are more likely to be in the definition of the ideal structural properties of the formed materials. Production of nitrogen is a kinetic separation that requires materials with small pore openings. Again, although this process is well-established, the data accumulated over the years can provide a benchmark to understand whether accurate a-priori predictions based on force fields that are efficient in equilibrium calculations can be used also in predicting diffusivities.
Finally, we envision that other separation processes, such as membrane separations, where the performance of the process is defined by the material used, will also benefit from multiscale screening workflows in producing more realistic, performance-based rankings of the available materials.
The role of ML methods will grow: It is already evident that the scope of multiscale screening methods will be expanding along with the range of available materials. This, combined with a large number of parameters, leads to multidimensional material-process configuration-performance space, which is more efficiently navigated using ML methods. This direction is both very promising and still widely unchartered. Hence, there is a strong incentive to more fully explore potential of the ML models to accelerate process-level screening workflows.
Quality data, reproducibility of results and consistency of comparison: we believe these aspects will be a singular, most important barrier for the multiscale approaches to make an actual impact through identifying both better and realistic options for carbon capture. The molecular simulation community has already produced a substantial number of screening studies for carbon capture. Similarly process level community has been examining various options for both processes and materials (but not in a large scale screening mode) for this task. The multiscale methods emerging from the combination of these two realms have been reviewed here. However, these studies use different assumptions, models and conditions which makes systematic comparison of their results difficult. One possible proposal for the simulation community would be an open call for the systematic comparison of the currently existing process modelling codes (including commercial ones) and model assumptions using a reference case study. This will be a significant step towards building confidence in ranking of the materials.
Techno-economic analysis: Development of multiscale screening studies should eventually go beyond the process-level. This is because the ultimate driver for commercialization of adsorption-based carbon capture processes is the cost. Therefore, the screening studies at the process level must be linked with techno-economic analyses where the ultimate design objective is to reduce the overall cost of CO2 capture and concentration (CCC). Although, there has been some recent attempts at this direction [52][53][54], there is still a dire need for integrated adsorbent-process optimization linked with techno-economic assessments of the CCC technology to justify its commercialization.
The ultimate challenge in post-combustion carbon capture still remains: It is important to recognize, that despite 15 years of computational material screening studies for carbon capture, an actual pilot scale plant based on a MOF or ZIF has not been built. While the aim of developing an in-silico route to finding optimal materials is a sound aspiration, there is the need to include in the selection process also the ability to synthesize the new materials and assess their stability against thermal cycling and exposure to contaminants and moisture. Predicting the stability of the materials is a challenging area. There are also other technical issues associated with scale-up of the process including the number of adsorption columns required and their control in large commercial plants, the vacuum level needed for a VSA process, the plant footprint as a whole, and finally the cost. Hence, it is clear that there are still significant challenges towards industrial implementation of carbon capture technologies based on novel materials. Notwithstanding, we believe development of more advanced mutliscale methods (such as those reviewed in this article) is an important step for accelerating our progress towards this objective.