Transferable force-field for modelling of CO2, N2, O2 and Ar in all silica and Na+ exchanged zeolites

In this work we propose a new force field for modelling of adsorption of CO2, N2, O2 and Ar in all silica and Na+ exchanged Si–Al zeolites. The force field has a standard molecular-mechanical functional form with electrostatic and Lennard-Jones interactions satisfying Lorentz–Berthelot mixing rules and thus has a potential for further extension in terms of new molecular types. The parameters for the zeolite framework atom types are optimized by an iterative procedure minimizing the difference with experimental adsorption data for a number of different zeolite structures and Si:Al ratios. The new force field shows a good agreement with available experimental data including those not used in the optimization procedure, and which also shows a reasonable transferability within different zeolite topologies. We suggest a potential usage in screening of different zeolite structures for carbon capture and storage process, and more generally, for separation of other gases.


Introduction
A rapid increase of greenhouse gases, particularly carbon dioxide, in the atmosphere is one of the actual global problems. Thus, a very important issue in modern society is to prevent CO 2 emission from different kind of sources into the atmosphere. Many efforts have been addressed particularly in the development of cost-efficient technologies of separation and capture of carbon dioxide [1] from the fuel gases. A general opinion is that Carbon Capture and Storage (CCS) process with a good choice of adsorbents can satisfy low cost energy requirements. [1][2][3] There is a large number of different materials which can be used as CO 2 adsorbents, and finding the most efficient adsorbents has attracted a lot of both experimental and theoretical research [4][5][6][7][8][9][10]. Zeolites is a perspective class of materials suitable for CO 2 capture due to their high adsorption capacity and easy regeneration.
Experimentally it might be difficult if not impossible at all to test many thousand of different zeolite materials and their variations which can be used as CO 2 adsorbents. Computationally based studies of CO 2 adsorption using molecular simulations can provide important information about how good are specific materials for CO 2 capture, and thus be used as a screening tool for the choice of the most efficient adsorbents. These studies however require an accurate force field which should have a good transferability in between different sorbent frameworks, in order to provide reliable and unbiased description of their adsorption.
Currently one can find in the literature several force fields suitable for studies of gas adsorption properties. Some of these force fields are of universal type, for example UFF [11], Dreiding [12], COMPASS [13]. These force fields include parameters for many different atom types, but they were not specially optimized for modelling of gas adsorption. Using them for modelling of experimental CO 2 adsorption isotherms produces often poor results. For the purpose of screening, it is especially important to optimize force field parameters for a certain class of material (e.g. zeolites) as well as for the target properties (adsorption load, adsorption enthalpy) in the relevant range of thermodynamic conditions. It is also necessary that the force field be transferable in between different types of zeolites, that is the force field with same set of parameters is able to reproduce experimental data for different zeolite topologies, as well as for varying composition, and not be biased to a limited set of structures used for the force field parametrization.
There are several force fields which were parametrized especially to match exper imental data in terms of uptake or self diffusion of adsorbed molecules in zeolite-like sorbents [14][15][16][17][18][19]. Almost all of those force-fields are derived for all-silica or alumino-silicate zeolites with sodium ions as cations for charge balancing. The number of considered zeolite structures was however rather limited, and parameters were often optimized for CO 2 adsorption only [17,20]. Force field parameters for other molecules exist as separate individual studies, and they are not 'merged' together into a single force field. Any decent try of this would lead to a significant reparametrization of the force field parameters, especially those which concerns Lennard-Jones interactions. In the force field of Jarmillo and Chandross [14] available molecules are CO 2 , NH 3 and H 2 O, and parameters for zeolite framework are optimized for pre taken parameters of quest molecules. Two water models are considered, SPC/E and TIP3P. This force field is optimized for adsorption in 4A zeolite (LTA zeolite with sodium cations), however, it was not tested for other frameworks. In another force field [21] two guest molecules are considered, CO 2 and water, where parameters for CO 2 -zeolite interaction where taken from García-Sánchez et al [17] and cross term parameters for water-zeolite interactions where obtained by a combination of the SPC/E water, CO 2 and framework parameters. Yet, by testing all the mentioned force fields against available experimental adsorption data, we found that non of these force fields are able to match available experimental data for large number of zeolites and several types of guest molecules. Moreover, most of these force fields introduce cross term interactions as independent parameters (without using Lorentz-Berthelot or other combination rules [22][23][24]) which leads to unnecessary large number of parameters and possible overdetermination of the optimization procedure, creating also difficulties in possible further extension of the force field.
Thus, an important and still open issue in computational zeolite science is to create a force field for accurate and predictive reproducing of experimental adsorption isotherms of different adsorbate molecules. Our work is motivated by this matter, and here we propose a new force field for adsorption of small molecules in Na + exchanged zeolites. Using this force field we are able to reproduce experimental data for adsorption of carbon dioxide, nitrogen, oxygen and argon in an representative set of all silica and Na + exchanged zeolites for thermodynamic conditions relevant to CO 2 separation processes. The force field covers a wide range of experimental adsorption data which we found in the literature. The force field was validated by computation of adsorption data which were not included into the optimization procedure. According to our results, this force field have shown an overall good transferability in between different zeolites, as well as within same zeolite topology but with varying Si:Al ratios. An important feature of the developed force field is that all cross-interactions between different atoms are determined using Lorentz-Berthelot mixing rules which simplifies possible extensions of this force field.

Force field description
The total energy in this force field is described in a standard way [12,13,17,25,26] as a sum of bonded and non-bonded pairwise interactions: whereas the bonded interaction is given as a sum of bond stretching and angle bending A pair of atoms is considered as non-bonded if they belongs to different molecules or if they are on the same molecule but separated by more than two covalent bonds. A case when two atoms are in the same structure and separated by three covalent bonds, so called 1-4 interaction, is treated with scaled Lennard-Jones and Coulombic potentials with scaling factor of 0.5. All 1-2 (directly bonded atoms) and 1-3 (two atoms separated by two bonds) non-bonded interactions are ignored, since they are included in the bonded part through bond stretching and angle bending. Since zeolites are crystalline materials with small fluctuations of the framework, dihedral potentials are not included in this force field. For zeolite framework, Dreiding force field [12] was used to describe bonded interactions.
For modelling of adsorption properties the most important part of the force field is the nonbonded interaction, which is given as a sum of steric and dispersive van-der-Waals interaction modelled by Lennard-Jones 12-6 potential and electrostatic interactions where σ ij is a distance at which the potential energy is equal to zero, ε ij is a depth of the potential curve, ε 0 is the vacuum permittivity ( × − 8.85 10 12 C 2 J −1 m −1 ), r ij is distance between sites, q i and q j are partial atomic charges of the given atom pair i, j.
For cross term pairwise Lennard-Jones interactions, Lorentz-Berthelot combination rules are used: Parameters of the non-bonded interactions (σ, ε and q) for each atom type were optimized in this work to provide the best possible matching with experimental adsorption data.
Zeolite structures used in simulations were prepared in a standard way [14,15,30]. The initial framework structures were taken from the IZA zeolite data base [31]. As a rule, all T-atoms in these structures were silicon. The unit cell was periodically repeated along X, Y and Z axes so that size of the supercell in each direction was not less than 24 Å (double of the cut-off distance equal to 12 Å). Then a certain number of silicon atoms within the zeolite super-cell were replaced by aluminum to have desirable silicon/aluminum ratio, which is the number of Si atoms divided by the number of Al atoms. Since positions of Al 3+ ions in a zeolite framework are unknown for most of the zeolite topologies, in this work we used a random method for distribution of Al 3+ ions under the only constrain that Löwenstein's rule [32] (preventing two Al atoms to be connected by a single O atom) should be satisfied.
Each time silicon was replaced by aluminum ion, zeolite framework became negatively charged. In order to balance the framework charges, sodium cations Na + were added into the cell in a number equal to the number of substituted Al atoms. A simulated annealing procedure with an energy minimizing algorithm [33] was used for this purpose. After that, zeolite lattice atoms were allowed to relax by energy minimization within the Dreiding force field [12] keeping the simulation cell geometry intact, in order to optimize the structure of zeolite atoms taking into account Si to Al substitution and inserted sodium ions. Typically we observed only a minor adjustment of the atomic positions (compared to the structure imported from the Zeolite database) as a result of the energy minimization procedure.
In this work we distinguish two types of zeolite oxygen atoms. The first one is oxygen atom attached to two silicon atoms, which is denoted as O_z, and the second oxygen type which is bridging Al 3+ and Si 4+ cations, denoted as O_z2 force field type, see figure 1. Both partial charges and Lennard-Jones parameters for these two types of oxygens were optimized separately in the optimization procedure.

Models for guest molecules
For guest gas molecules we used force field parameters available from the literature. We have chosen models (for CO 2 , N 2 , O 2 and Ar), figure 1, which were proven to reproduce bulk properties of gases, their liquid state and gas-liquid equilibria. We have controlled the considered models by reproducing liquid-vapor coexistence curve and we present results of these tests in the supplementary information, figures S1-S3 (stacks.iop.org/MSMSE/24/045002/mmedia).
Carbon dioxide molecule was modelled as a three site molecule with parameters from papers [34,35] (EPM2 model). This model was parametrized to reproduce properties of the gas-liquid coexistence including the critical point region. Carbon atom in the middle of the molecule (C_co2) is connected to two oxygen atoms (O_co2) by chemical bonds 1.149 Å long. Although CO 2 molecule is a small linear molecule and thus often modelled as a rigid one, in this work it was modelled as flexible with parameters taken from the original paper [34]. We have also tested a few other CO 2 models, particularly rigid EPM2, EP [34] (a model preceding EPM2) and TraPPE [36], and we provide results of this testing in the supplementary information, figures S1, S4 and S5.
Nitrogen molecule N 2 was modelled as a three sites model [37] which consists of two positively charged nitrogen atoms (N_n2) with inter atomic distance of 1.098 Å, connected via massless virtual site N_mid. A positive massless charge point q + , equal to | | − q 2 (where q − is a partial charge of nitrogen atom), was placed at the molecular centre of mass. In this way the experimental quadrupole moment of an N 2 molecule was emulated, which is an important physical property [1,2,38] affecting gas adsorption [39]. Representing N 2 as a two site model with zero partial charges leads to strong underestimation of N 2 adsorption in comparison with experimental data even if liquid-vapour coexistence curve is reproduced well (figures S2, S6, S7). The Lennard-Jones parameters of N 2 model were tuned to reproduce liquid N 2 properties by Murthy et al [37], and the model have been used in a number of adsorption studies [40,41]. N 2 molecules were considered to be rigid in our simulations.
Although reasonable good agreement with experiment for the adsorption isotherms with the original force field parameters for N 2 was obtained, we have also considered further optimization of Lennard-Jones parameters for N 2 which showed a better agreement with available experimental adsorption isotherms. The details of re-optimized N 2 parameters and compariso n between adsorption isotherms from experimental data and computed with re-optimized forcefield are discussed further in the text.
Oxygen molecule O 2 was modelled similarly to nitrogen as three sites molecule with two oxygen atoms (O_o2) connected via massless point charge (O_mid). Lennard-Jones parameters and partial atomic charges for O 2 molecule were taken from the works [42,43]. Presented 'molecules' are: zeolite 4 member ring, sodium ion, carbon dioxide, nitrogen, oxygen and argon. Codes denote force field type assigned to particular atoms.
Inter-atomic distance between two oxygen atoms was 1.017 Å. During the simulations, oxygen molecules were treated as rigid.
Argon is modelled as a single spherical atom site with zero charge and Lennard-Jones parameters taken from [44].

Computation of adsorption isotherms
All adsorption isotherms in this work were computed by the Grand Canonical Monte Carlo (GCMC) method for gas molecules in a supercell of zeolite sorbent with periodic boundary conditions. Within GCMC, by random insertion and/or deletion of guest molecules with probability of acceptance determined by the rules of the Grand Canonical ensemble, average number of molecules inside the zeolite framework can be determined for a specified value of the chemical potential. CGMC simulations reproduce thus the average adsorption load of the gas molecules in the bulk zeolite material which is in thermodynamic equilibrium with gas molecules having the same value of chemical potential which in turn is defined by the gas fugacity. In this work we used the Peng-Robinson equation of state [45] for conversion between gas pressure and fugacity. Most of simulations have been carried out by 'RASPA' software [46][47][48]. We have also used software packages 'Music' [49,50] and Materials Studio Accelrys [51,52], essentially the same results (under the same settings) were obtained with all three packages.
In all simulations the total number of MC cycles was taken to be at least × 5 10 5 , where initial one third of total steps was taken for equilibration and the rest of the simulation for collection of averages. The Lennard-Jones interactions between particles were calculated within cut-off radius equal to 12 Å, which was less than a half of the supercell unit cell vector lengths in all our simulations. The potential was truncated at the cut-off distance and long range tail correction was used. In principle, use of tail corrections allows for even shorter values of the cutoff for zeolites with small unit cell (which was used in some previous simulations), but it is worth to note that EPM2 CO 2 model was optimized with the same settings [34] as used in the present work. The Ewald summation method [23,24,53] was used for treatment of the electrostatic part of the interaction potential.
Many zeolite structures contain small cages (sodalites) surrounded by 4-or 6-rings which are not penetrable by gas molecules. These cages need to be blocked in order to prevent possible insertion of a guest molecule inside such cages within the Grand Canonical procedure. We block inaccessible pockets in zeolite framework utilizing our code designed for this purpose and validating results with those obtained using Zeo++ software [54][55][56]. Further details of blocking procedure as well as investigation the effect of blocking within the developed force field is given in the supplementary information, figures S8-S11.
During grand-canonical simulations of the adsorption isotherms, the framework atoms and blocking pseudo atoms were kept fixed at their initial crystallographic positions, while sodium ions are allowed to move. Analysis of the sodium 'mobility' reveals that sodium ions are typically located near positions determined in the optimization procedure with deviations within 1 Å. To compare results from simulations with fixed and free sodium ions, we run a short molecular dynamics (MD) simulation with 2 fs time step and 10 5 MD steps in NVT ensemble and collected 10 random configurations from the MD trajectory. For MD simulations we used MDynaMix Package [57]. GCMC simulations were then performed on all 10 samples with fixed ions and average value of them was compared with the GCMC data obtained where Na + ions were allowed to move. The difference was found not to exceed 1-2 percent being of the order of the statistical uncertainty (see figures S12 and S13 of the supplementary information).
We have also tested the role of flexibility of the framework and of the CO 2 sorbent molecules on the adsorption results. It was previously reported [30] that in case of bigger sorbents (such as octane, benzene, etc) flexibility of the zeolite framework may be of importance, while for smaller molecules inclusion of the framework flexibility does not affect the adsorption results. Our calculations shown in figures S14 and S15 of the supplementary information demonstrate that CO 2 adsorption curves for FAU (Si:Al = 1.4) and LTA (Si:Al = 3.5) zeolites computed within rigid and flexible framework models (the later described by the Dreiding force field) coincide within less than 1%. Comparison of CO 2 adsorption in FAU (Si:Al = 1.4) framework computed with rigid and flexible CO 2 models shows that flexibility of the CO 2 model does not affect the adsorption results either. More details on the role of flexibility and other details of the molecular models involved in the simulations are given in supplementary information, sections 2-7.

Optimization of the force field
In developing of a physically grounded force field, one usually relays on results of quantumchemical calculations or on experimental data (or on combination of both). A major obstacle in parametrization of a force field for adsorption in zeolite-like structures is that quantumchemical methods are hardly applicable, mainly because of importance of van-der-Waals (dispersion) interactions which are generally poorly described within DFT approaches, uncertainties in determination of partial atom charges in periodic systems and other computational limitations. Therefore in this work we used available experimental data as the main source of information for the parametrization of the force field. Similar approach have been used previously in literature [14,17], where the force field parameters were fitted to adsorption isotherms data measured in experiments. There are however two major differences from previous works. First, here we perform an optimization of the force field parameters to a large set of adsorption data which includes different framework types, different Si/Al ratios and different temperatures. Second, we use the same set of parameters for all zeolite frameworks to describe adsorption of different molecules, using Lorentz-Berthelot combination rules for cross-interactions.
The framework structures for which computation of the adsorption isotherms have been carried out, with references to the works with experimental adsorption data for the same structures, are listed in table 1. This list is divided in two parts. The first part, containing 11 structures (28 computations in total accounting for different temperatures), was used for the optimization of the force field parameters. The second part, including 10 structures, was used for the validation of the derived force field.
An automatized iterative approach was used for tuning the parameters of zeolite atoms based on minimization of the relative error in comparison to experimental data of pure CO 2 , while parameters for CO 2 atoms were taken from literature [34,35]. First, we made an initial guess for the force field parameters for zeolite framework. This was done by taking approximate average values of the Lennard-Jones parameters and partial atomic charges from the literature [14][15][16][17][18]21] with exception for sigma parameters of Al and Si atoms, which are taken from Dreiding force field [12]. Next, we selected the set of experimental adsorption isotherms (part A of table 1) and on the basis of the initial force field guess we performed GCMC computations at least in five fixed pressure points for each of the isotherms. If an experimental adsorption isotherm showed significant uptake at initial pressure stage, we used a logarithmic scale for the choice of the five equally distributed pressure points in which we compared the computed and the experimental data, otherwise (with a Henry law-like isotherms) we use a linear pressure scale. We compared computed adsorption loading with the corresp onding experimental data and calculated an average error as a difference between computed and experimental data, by using the following expression exp. (7) where N is the total number of fixed pressure points computations for all isotherms in a single iteration, n exp. and n sim. are corresponding experimental and simulated loadings, respectively. The convergence criteria was (somewhat arbitrary) set to 5% of the average deviation between corresponding experimental and computed values of the adsorption load. If the calculated average error exceeded the predefined convergence criteria, an update of the force field parameters was carried out. While optimizing parameters, care was taken to have a set of the force field parameters with reasonable physical meaning which were also consistent across the whole range of Si:Al ratios. Particularly we were restrictive to change Lennard-Jones parameter sigma expressing effective atom diameter. We have used fixed charge 1.0e for Na + counterions (some previous force field [17,[68][69][70] used scaled values of Na + charges referring to the charge transfer effect). Furthermore, while varying partial charges of the framework atoms, overall electroneutrality should be kept which put some constraints on possible values of the partial charges. For all-silica structures with only two types of framework atoms this constraint is evident: q(Si_z) = −2q(O_z). Each time when an Al atom substitutes a Si atom, additional Na + counterion is inserted, and four oxygen atoms change their type from O_z to O_z2. To satisfy the electroneutrality condition at Si to Al substitution, another constraint should be put on the partial charges: q(Al_z) + q(Na +) + 4q(O_z2) = q(Si_z) + 4q(O_z). These constraints  [67] were satisfied in the optimization procedure, and they also hold for the finally optimized values in table 2, providing thus consistent description of electrostatic interactions in zeolite frameworks of any Si:Al ratio. Each time when the program have to update the force field within the iterate procedure, test computations were carried out with varied values of each of the parameters in order to understand how each parameter affect the adsorption loading of pure gas component inside a zeolite framework. An example of such series of computations carried out with initial force field parameters and scaling of either charges or parameter Lennard-Jones epsilon is shown in figure 2. One can see that increase of both particle charges and Lennard-Jones parameters ε lead to increase of the adsorption load, but the shape of the adsorption isotherm changes differently with increase of either q or ε.
After each iteration the difference in between computed and experimental values at high and low pressures is evaluated, and a heuristic script (more advanced optimization procedure can be found in papers [71,72]) suggests changes of parameters which would provide better agreement between computed and experimental adsorption isotherms. On the initial stages of the iterative procedure a general scaling of either partial charges or Lennard-Jones parameters were done, in the later stage tuning of parameters of individual atom types (Si, Al, and oxygens of the framework) was done to satisfy individual isotherms of the considered set of structures. The finally optimized values of parameters are listed in table 2. Note: The first column is a molecular type for which force field is given. The second column gives notation of the force field atom type. The description of the force field atom type is given in the third column. Fourth, fifth and sixth columns of the force field table stand for sigma (Å), epsilon (K ) and partial atomic charge (e) for the given atom type, respectively. Seventh column is the reference to the work from which parameters were taken.
In our computations we found that shape of the adsorption isotherm curves is strongly influenced by the Lennard-Jones parameters of aluminum and silicon atoms. These parameters were not included into many previous force fields, with argumentation that since the polarizabilities of aluminum and silicon are lower than those of oxygen atoms, it would make a decent assumption to ignore Al and Si atoms in Lennard-Jones interactions and thus include only oxygen, as in the work of García-Sánchez et al [17], Newsome et al [68] and Maurin et al [73]. However, this assumption in many cases tends to underestimate adsorption at low pressures and also overestimate adsorption of CO 2 at high pressures [17,29]. Explicit consideration of Lennard-Jones interactions at T-atoms gives better results in terms of transferability in between different zeolites, thus they were included into our force field.
Using this methodology we were able to derive the force field which is able to reproduce experimental data for many different zeolite structures and thus obtain a good transferability of the force field. The predictive character of the derived force field was confirmed by a series of computations with zeolite structures not included into the optimization procedure (part B of table 1). As it will be shown later in the validation section, using the optimized parameters for zeolite framework (obtained by optimizing the interaction with CO 2 only) we were able to reproduce most of the experimental data for N 2 , O 2 and Ar adsorption. This gives an additional argument to the validity of the present force field.

CO 2 adsorption in all-silica zeolites
We begin comparison of simulated and experimental adsorption isotherms of pure carbon dioxide for all-silica structures. In the first example of all silica CHA zeolite, we compare our simulation data with experimental data taken from work of Pham et al [41], see figure 3(a). The adsorption isotherms are in this case measured in a wider range of temperatures from 273 up to a 343 K. In this example our force field tends to underestimate adsorption loading at higher temperatures. Note however some differences between experimental data of works [61] and [41] for the same CHA structure at 301 and 303K.
We compare also our simulated data with the experimental data from work [41] for all silica FER zeolite, see figure 3(b). Adsorption isotherms are measured and computed at 273, 303 and 343 K. As in the previous case, in this example simulated data are somewhat underestimate the experimental adsorption loadings.
In the next example we compare our data with adsorption isotherms for all silica BEA zeolite taken from the same work of Pham et al [41], see figure 4(a). Adsorption isotherms are measured at 273, 303 and 343 K. Unlike to previous examples, in this case our force field was able to reproduce these experimental data only qualitatively, with underestimation of the adsorption loading. We include also in figure 4(a) adsorption isotherms data computed with previous force field from the work of Makrodimitris et al [35]. Still the force field developed in our work demonstrates certainly an improvement in comparison to the previous force field.
In figure 4(b) we displayed the isosteric heat of adsorption of CO 2 in zeolites measured in paper by Pham et al [41] at 303 K in comparison with simulated data. A good agreement is seen for all silica CHA, BEA and MFI zeolites, while for CHA zeolite an opposite trend in the behaviour of adsorption heat as a function of adsorption load is observed. The experimental trend for FER structure is also opposite to other measured structures, and may be caused by some reasons specific for FER structure.
To conclude with all silica zeolites, we compare simulated adsorption isotherms with experimental data from Sun et al [62] ( figure 5(a)) and García-Sánchez et al [17] ( figure 5(b)) for all silica MFI zeolite. In the work of Sun et al, the adsorption loading of pure CO 2 was measured in a wide temperature range from 276.6 up to 352.6 K, while the pressure range was up to 2000 kPa. Simulated adsorption isotherms are in a very good agreement with experimental data in the whole pressure and temperature range. The same all-silica MFI structure was studied experimentally in the work of García-Sánchez et al [17]. Comparison of our results with that work is shown in figure 5(b). Also shown are results of computed adsorption isotherms using the force field developed in paper [17] figure 5(a)) where used in the validation while experimental data from García-Sánchez et al [17] ( figure 5(b)) where used as a part of the force field optimization procedure.

CO 2 adsorption in zeolites with varying Si/Al ratio
In another series of examples, we compare simulated and experimental data for zeolites with the same structure code but varying Si/Al ratio. In figure 6 we compare computed adsorption  isotherms in FAU zeolite with three different Si:Al ratios, with experimental data taken from papers Dunne et al [64] for Si:Al = 1.23, Akhtar and Bergström [58] for Si:Al = 1.4 and García-Sánchez et al [17] for Si:Al = 2.5. Adsorption temperatures (273, 303 and 304 K) are indicated at the corresponding adsorption isotherms in figure 6. Our computed data shows very good matching with experimental data from papers [58] and [64], while in comparison to data from [17] the computed adsorption isotherm somewhat overestimates the experimental  [41], previous force field (lime coloured closed symbols) [35] and this force field simulation data of adsorption isotherms of pure CO 2 in all silica BEA zeolite. Adsorption isotherms are measure and computed at 273, 303 and 343 K. (b) Comparison of experimental [41], and computed isosteric heat of adsorption at 303 K, for all silica CHA, FER, MFI and BEA zeolites. 273 K -sim. [38] 303 K -exp. 303 K -sim.
343 K -sim. [38] (a) adsorption loading. Anyway simulations reproduce well the trend of increase of CO 2 adsorption with decrease of Si/Al ratio in the FAU zeolite structure.
The adsorption uptake of pure CO 2 in MFI zeolite with Si:Al ratio 130 was measured at temperatures 308, 323, 358 and 393 K in the work of Ohlin et al [59]. In comparison to this set of data we found that at adsorption temperature 308 K the simulated isotherms are in a very Experimental data from García-Sánchez et al [17] (open symbols), data from previous force field [17] (lime colored symbols) and simulated isotherms with this force field, for temperatures 253, 273 and 303 K.   253 K -sim. [17] 273 K -exp. 273 K -sim.
303 K -sim. [17] (b) good agreement with experimental data, figure 7. However, at higher temperatures our force field overestimate the experimental adsorption loadings. The maximal discrepancy in between simulated and experimental adsorption loadings can be obtained at temperature and pressure equal to 358 K and 100 kPa, respectively.

Validation simulations
A real test of a newly derived force field is to match the results to experimental data not used for the optimization of the force field parameters. In figure 8(a) simulated adsorption loading is compared with experimental data from work of Fang et al [61] for all silica CHA zeolite, measured at 301, 313 and 323 K. Fang et al in the same work derived force field parameters describing CO 2 adsorption from periodic dispersion-corrected DFT calculations for all silica structures, and carried out validation computations only for all silica CHA zeolite. They obtain perfect matching with their experimental data. Our force field also shows almost perfect matching figure 8 except slight underestimation of adsorption at lowest temperature 301 K. Taking into account that our force field is not developed specifically for this zeolite structure, but rather for all zeolites with varying Si/Al ratio, the agreement with experimental data of work [61] is very good.
In order to validate representativity of the new force field (that is ability of the force field to reproduce other properties than adsorption loading), we compared, for the same all-silica CHA zeolite, the computed isosteric heat of adsorption with experimental data from the same work of Fang et al [61]. Our data for the heat of adsorption at 301 K shows a good agreement in comparison to experimental measurements, figure 8(b), although at high loadings, computed adsorption enthalpies are slightly lower that those in experimental measurements. According to works [2,[73][74][75][76], a probable explanation of this is slightly underestimated Figure 6. Adsorption isotherms of CO 2 in two NaX and one NaY zeolite structures (FAU structure code) with different silicon over aluminum ratios, Si:Al = 1.23 and Si:Al = 1.40 and Si:Al = 2.56. Experimental data for comparison are taken from [64,58] and [17], for ratios 1.23, 1.40 and 2.56, respectively. Open symbols refers to experimental data and closed symbols to simulated adsorption isotherm. adsorbate-adsorbate interactions at higher loadings (in this case the interaction between adsorbed CO 2 molecules), in comparison to experimental data.
In figure 9 we compare CO 2 adsorption data with experimental results for LTA zeolites of work by Martin-Calvo et al [20] (Si:Al = 1) and Palomino et al [60] (Si:Al = 3.5 and all silica), and for MFI zeolite with Si:Al = 280 from data by Harlick and Tezel [67], see figure 9. Agreement for all-silica LTA zeolite is perfect. For two other structures simulations somewhat underestimate adsorption, though overall agreement can be characterized as good. Comparison of simulated and experimental result for LTA Si:Al = 1 zeolite from work [20] is also presented in logarithmic scale in the supporting information, figure S16.
Further comparisons were done with experimental data for the MFI zeolite with Si:Al ratio 30 [64] and 120 [66], as well as for CHA structure with Si:Al ratio 6 from Pham et al [63], figure 10. Again, we see very good overall agreement, with almost perfect reproduction of experimental data for CHA structure. For MFI structure with Si:Al ratio 30, simulations somewhat underestimate adsorption at low partial pressures, while for MFI structure with Si/Al ratio 120 simulations slightly overestimate adsorption at higher partial pressures of CO 2 . Another comparison for MFI zeolite with Si:Al ratio 13 is shown in supplementary information, were adsorption isotherms computed with our force field are shown together with computed and experimental results from paper of Newsome et al [68].

Adsorption of other gas molecules
Following the procedure for validation of the force field parameters for adsorption of carbon dioxide, we tested the developed force field for other gases: nitrogen (N 2 ), oxygen (O 2 ) and argon (Ar). However, in comparison to available experimental data for adsorption of carbon dioxide, there is less data for adsorption of other gases at room temperatures, where the major part contributes from nitrogen. Considering separation of gases using zeolites as adsorbents [1,3,[77][78][79] or membrane based separation [80,81] one can state that knowledge of  interactions between other than CO 2 guest molecules and zeolite framework is of importance for this class of applications. Nitrogen adsorption. Using partial atom charges (reproducing the N 2 quadrupole moment) and Lennard-Jones parameters for nitrogen molecule from the literature [37], as described above, we where able to reproduce quantitatively most of the experimental data availavle  for N 2 adsorption at ambient temperatures. In the first example of validating force field parameters for nitrogen adsorption, figure 11(a), we compare experimentally measured adsorption isotherms of pure nitrogen, at 303 K, for all silica FER, CHA, MFI and STT zeolites, with the corresponding simulated data. The agreement with experiment can be judged as very good, with simulation results slightly underestimating experimental adsorption for MFI and STT zeolites.

(b)
In the next example we compare computed and experimental adsorption isotherms for pure nitrogen at several temperatures in two zeolite structures with high Al content, figure 11(b). While overall behaviour of the uptake with respect to N 2 pressure is found to be good, the simulation results underestimate adsorption by factor about 1.7.
Comparison of simulated and experimental adsorption isotherms measured at different temperatures for FAU zeolite with Si:Al = 1 from paper by Park et al [83] is shown in figure 12. The agreement is better than in the previous case, with underestimation of the simulated adsorption load by 15-20% of the magnitude.
Taking into account that the amount of the adsorbed nitrogen in all examples is low (in comparison with typical CO 2 adsorption), the overall agreement between simulated and experimental adsorption data can be considered as satisfactory. Still, for more reliable prediction of the competitive CO 2 -N 2 adsorption for screening purposes, additional optimization of nitrogen force field parameters may be desirable. In the example above we used N 2 parameters optimized for properties of liquid nitrogen [37] and they might not be optimal for N 2 gas phase at temperatures and pressures typical for conditions of CO 2 separation process. We therefore reoptimized parameters for nitrogen molecule, following essentially the same procedure as for optimization of zeolite framework parameters. In this procedure, only LJ parameter ε was varied in order to fit adsorption isotherms of works [41,64,82,83]. Value / ε = k 40.24 b K found to reproduce the experimental data in the best possible way as shown in figures S18-S20 of the supplementary information. It still should be noted that the reoptimized nitrogen parameters are verified in case when adsorption temperature is around 300 K and adsorption pressure is below 1 bar, and they should be used with care for modelling of adsorption at lower temperatures and/or higher pressure.
Oxygen adsorption. In comparison to previously reported data for carbon dioxide and nitrogen gases, there is significant lack of experimental data for oxygen gas. Thus, we validate our force field with the set of available experimental data. Knowing the force field parameters   for oxygen-zeolite interaction is important for several reasons. First of all, oxygen is quite often a part of the flue gases [7,84] thus knowing how this molecule interacts with zeolite framework is of great importance for modelling of gas separation using zeolites as adsorbents. Furthermore, since oxygen is quite often separated from atmosphere components [85] and used in medical purposes, without a good force field for adsorption of this gas, it is not possible to screen different types of zeolites for separation of oxygen from other components. To the best of our knowledge, we are not aware of such a study, so it would be interesting to see how our force field perform in that type of application.  [64] for Si:Al = 1.23, Jayaraman et al [82] for Si:Al = 1.25, and Park et al [83] for Si:Al = 1.
Comparing experimental adsorption isotherms with computed in this work using new force field parameters for zeolite framework atoms and oxygen parameters from the literature, we found a very good agreement with available experimental data figure 13. This validation, together with shown in the previous section results for nitrogen adsorption, gives a strong argument for the physical consistency of the present force field.

(b)
Argon adsorption. Similarly to oxygen, the presence of argon quite often can be seen in the flue gases, as well as a minor part of the atmosphere. This gas, however, is not so attractive in separation from other components as a primary component, as CO 2 might be. Instead, the presence of argon affects separation of other gases from the gas mixture. Adsorption of argon is usually measured in experimental studies for evaluation of the surface area, which is done at low temperatures. This low temperature adsorption isotherms may be however not well suitable for validation of the force field parameters for zeolite-argon interaction at ambient temperatures typical for fuel gases. We therefore used a few available literature data for Ar adsorption in zeolites at room temperatures. Comparison of computed and experimental adsorption isotherms of pure argon gas in FAU zeolite with several different Si:Al ratios as well as in MFI all silica zeolite is shown in figure 14. Experimental data for adsorption of argon in NaX zeolite are reproduced almost perfectly, while computed adsorption isotherms for MFI zeolite overestimate slightly the experimental data.

Modeling of gas diffusion
Finally, we demonstrate suitability of the developed force field for molecular dynamics simulations and studies of kinetic properties of the adsorbed molecules. We have carried out molecular dynamics simulations of CO 2 molecules adsorbed in all-silica MFI and in FAU (Si:Al = 1.23 ) framework structures (known also as ZSM-5 and NaX zeolites) and computed self-diffusion coefficients of the guest molecules. In the simulations, the unit MFI cell was periodically replicated × × 3 3 4 times and FAU cell × × 2 2 2 times, providing the final supercell of about 50 Å in each direction. The supercells were loaded by CO 2 molecules counting 7 molecules per 48 framework T-atoms, corresponding approximately to 2.1 and 2.4 g mmol −1 load for FAU and MFI respectively. The cell geometry was fixed while all the atoms including the framework were allowed to move. Simulations were performed by MDynaMix v.5.2.7 package [57] using the double time step algorithm with short time step of 0.2 fs for intramolecular (CO 2 ) and framework vibrations, and 2 fs for intermolecular  interactions. The electrostatic interactions were treated by the Ewald summation method. The starting configurations, obtained from MC simulations, were preequilibrated during 300 ps at 150 K, then 6.5 ns runs were performed at each temperature point. Initial 1.5 ns of the production simulations were omitted from the mean square displacement analysis and remaining 5 ns of the trajectories were used for self-diffusion coefficients computations from the slope of the linear part of the mean square displacement curve according to: The temperature dependence of the computed self-diffusion coefficients is displayed in figure 15 in a form of Arrhenius plot, together with experimental data from paper by Kärger et al [86]. One can see a perfect match for the CO 2 diffusion coefficients in FAU zeolite structure. For MFI zeolite simulations show somewhat faster diffusion than experiment, nevertheless the slope of the simulated and experimental dependencies is the same, meaning that the activation energy of barrier crossing is reproduced correctly. These results validate further the developed force field and demonstrate its applicability for studies of adsorption kinetics.

Conclusion
In this work we have developed a new force field for accurate description of the adsorption properties of several gas species in Na + exchanged zeolites. As we have shown, this force field gives a very good transferability between different zeolites, as well as within same zeolite framework type with different Si:Al ratios. The force field is derived in an iterative procedure where the initial force field parameters for zeolite framework atom types are mostly taken as an average value of those in the present literature, and optimized to fit the available adsorption data. During the optimization of the force field parameters, the difference between experimental and computed adsorption loadings for different molecules was systematically minimized. An important feature of the minimization is that the same set of force field parameters was fitted at once for a large set of experimental data including different framework types, Si/Al ratios and temperatures. This procedure provided a good transferability of the derived force field, which was demonstrated by validation simulations of several zeolite structures not included into the fitting procedure. The difference between experimental and computed adsorption is often of the same order or less than difference between experimental data obtained for the same zeolites but coming from different experimental groups.
We have also further validated the present force field by computations of adsorption of nitrogen, oxygen and argon gases in zeolites. The same parameters for the zeolite framework were used as for CO 2 adsorption and parameters for all guest molecules were taken from the literature, which gives a strong evidence of the physical consistency of the present force field. The finally optimized Lennard-Jones zeolite parameters and partial charges of the zeolite framework are similar to commonly used force fields for the considered atom types which gives confidence that the optimized parameters are realistic and can be used to simulate adsorption of gas mixtures and adsorption selectivity. We have also demonstrated suitability of the new force field for modeling of kinetic properties by molecular dynamics simulations.
The validated transferability of the newly parameterized force field and its ability to reproduce adsorption of several gases makes it suitable for screening of different zeolite structure for the most efficient sorbents for CO 2 capture, as well as for other gas separation applications.