Isotherms of Fluids in Native and Defective Zeolite and Alumino-Phosphate Crystals: Monte-Carlo Simulations with “On-the-Fly” ab initio Electrostatic Potential

We use periodic Density Functional Theory (DFT) method to generate the electrostatic potentials of adsorption materials and use them in Grand Canonical Monte Carlo simulations of fluid adsorption isotherms. This permits us to consider complex solids showing defects and without a priori knowledge of their electrostatic parameters for the Monte Carlo simulations. We apply the method to aluminophosphate and silicate solids of ZON type and evaluate their affinity to adsorb and separate CO2 -N2 (H2 O) mixtures.


INTRODUCTION
Molecular mechanics methods have shown their efficiency at complementing and further validating experimental adsorption data of fluids in microporous solids [1,2]. However, the determination of the forcefield parameters that are used to compute the guest-host potential energy can be a difficult and time consuming process, especially for the parameters describing the electrostatic interactions. To overcome these difficulties, the electrostatic potential of the solid into which the fluid adsorbs can directly be used from periodic ab initio calculation [3]. This operation is done in a single step and only requires the solid crystal atomic positions. Another huge benefit of this approach is its versatility: there are virtually no limitations related to the solid crystal structure and the presence of defects in it (apart its size, which should be around several hundred atoms). To illustrate the strength and elegance of this approach, we consider the case of silicate and aluminophosphate crystals and the inclusion of silica defects in them, compensated or not by Na + . The CO 2 , N 2 and H 2 O isotherms in the native and defective crystal structures are determined and compared. The case of the CO 2 -N 2 mixtures isotherms is then investigated. Aluminophosphates (ALPO) have been shown to be interesting materials to separate CO 2 from CO 2 -N 2 mixtures as obtained from flue gas [4]. CO 2 separated from the mixture can afterwards be sequestered in geological deposits [5].

Solid Models
The starting ZON coordinates were taken from Roux et al. [6]. ZON has networks of parallel and zig-zag 8-membered ring micropores that interconnect and a zig-zag 6-membered ring network that also interconnects the large parallel network (Fig. 1) [7]. The symmetry group of ZON (1 and 2) is PBCA (orthorhombic+D 2h axis). The all silica equivalent is obtained by substituting all Al and P with Si. The silica defect ALPO models are defined by substituting twice adjacent Al-P pair (Fig. 1c) and a P atom (Fig. 1d). We sampled the different possible sites of Al-P pair substitution with Si-Si and computed their energies using MOPAC [8]. We selected the most energetically stable site. For material 4 (Tab.1), the electroneutrality of the system is verified after addition of Na + near the Sidefect and in the parallel 8-membered ring micropore (Fig. 1). Because the substitution of an Al-P pair with Si-Si pair does not change the global charge of the solid, charge compensating anion or cation is not needed. The introduction of defect in the ZON periodic unit cell lowers its symmetry group and the symmetry group is lowered to P1 for 3 and 4. Therefore, symmetry group P1 is used for all systems, including 1 and 2.

Periodic Density Functional Theory Calculations
The periodic DFT code VASP (Vienna Ab-initio Simulation Package) [9,10] is used to optimize the cell size and shape as well as all atomic positions with the only constraint that the cell symmetry remains orthorhombic, i.e., α = β = γ = 90º. The geometry optimizations are done using the conjugate gradient method as implemented in VASP-MedeA (Materials Exploration and Design Analysis) v. 4.6. The criteria to end the geometry optimization are considered met when all forces acting on atoms are less than 0.0193 kJ.mol -1 .pm -1 (0.02 eV.Å -1 ). The DFT functional that is used in the energy and force calculations is PBE [11,12] together with the Projected Augmented Wave (PAW) approach [13,14]. An energy cutoff of 38.596 MJ.mol -1 (400 eV) is set for the electronic wave functions. The optimized parameters of the cells are given in Table 1. After the cells are optimized, the total electrostatic potential calculations are done and the obtained electrostatic potential VASP files are converted to the format of the Monte Carlo code Gibbs without further modification.

Molecular Mechanic Monte-Carlo Simulations
All isotherm calculations are done using the Monte-Carlo code Gibbs v. 9.2 [15] as implemented in MedeA v. 2.8 [16].
Rigid molecule model is used for the fluid molecules, namely N 2 , CO 2 and H 2 O, in their pure or mixture state. For N 2 and CO 2 , the Lennard-Jones parameters, charges and geometries as given from the database of MedeA are used [17,18]. Water is modeled using the TIP3P model from Jorgensen et al. (1983) [18,20]. Its simulated saturation pressure is measured to be 27 mbar while the experimental value is 32 mbar at 298 K [19].  For the adsorption isotherm simulations of fluid on a solid, the solid is assumed rigid in the Gibbs code and Lennard-Jones and electrostatic potential energy grids characterizing the solid are generated prior the Monte-Carlo sampling tasks take place. The grids were generated for 2 × 2 × 2 extension of the primitive cells, i.e., the grid dimensions are ca. 2.8 × 3.0 × 3.4 nm 3 . In accord with previous results [20], only oxygen atom parameters are considered to generate the solid Lennard-Jones potential energy grid. The Lennard-Jones parameters for aluminophosphates [21,24] and silicates [22] oxygen atoms are different yet defined from the same group for consistency. They have been validated vs experimental data of alkanes and water adsorption isotherms [21][22][23][24][25][26][27].
All Monte-Carlo simulation runs are done in the Grand Canonical ensemble [24] at T = 298 K. For every simulation of a given fluid composition (pure or mixture) and fugacity, an equilibration run of five million steps takes place prior a production run of twenty five million steps. To ease the Monte-Carlo sampling, a configurational bias is used [25]. The Lennard-Jones pair energy cutoff is set to 2 nm, while the electrostatic potential energy is evaluated via Ewald summation, the optimal Ewald summation parameters being set automatically [30].
Prior to the determination the adsorption isotherms of N 2 , CO 2 , H 2 O and their mixtures, we determine the micropore volume to compute the solid porosity. This is done by helium-pycnometry simulations at low fugacity (0.001, 0.01 and 0.1 bar). At these low fugacities, helium behaves like an ideal gas and the density of helium in the solid is considered to be the same than that of bulk helium under the same conditions of fugacity and temperature. Hence, the determinations of the amount of helium adsorbed in the solid and of the bulk helium density give a direct estimate of the pore volume. Bulk helium density is determined by Grand Canonical Monte-Carlo simulation using the same temperature and fugacity than that used in the helium adsorption on the solids. The pore volume divided by the total solid volume gives the solid porosity. The porosities are reported in Table 1.
The water isotherm data are total number of adsorbed molecules, while that of N 2 and CO 2 are Gibbs excess adsorption values. The molecular excess adsorption is defined as: where n tot is the total number of molecules adsorbed on the solid at a given fugacity and composition in a given solid volume and n id , the number of molecule, non-interacting with the solid at the same fugacity and composition and in a volume equals to that of the micropore volume of the solid as determined by simulated helium-pycnometry. The n id values are determined using the same protocol than that used to determine the ideal helium density in the simulated helium-pycnometry. The excess adsorption values, N ex , are then expressed in g.cm -3 following: where M is the mole mass of the considered fluid (g.mol -1 ), V s the total volume of the solid for which n ex is measured (m 3 ) and N A the Avogadro number (mol -1 ).

RESULTS AND DISCUSSION
2.1 Influence of the Lennard-Jones and Electrostatic Potential on the Adsorption Isotherms of N 2 We consider the case of AlPO 4 and SiO 2 ZON frameworks before to simulate the adsorption of pure compound and mixtures of H 2 O, N 2 and CO 2 in the defective ZON solids. We limit the adsorption isotherm simulations to N 2 and analyze the consequences of changing the electrostatic and Lennard-Jones parameters. We calculate the Gibbs excess adsorption of N 2 between 1 and 120 bar and the isosteric heats of N 2 adsorption on the solids between 1 mbar and 1 bar using the AlPO 4 and SiO 2 frameworks, silicate and aluminophosphate Lennard-Jones oxygen parameters and electrostatic potential grids obtained from "classical" point charges [21,24,23,27] and ab initio calculations using VASP. The Lennard-Jones parameters for silicate and aluminophosphate oxygen atoms are 517.6 and 779.4 J.mol -1 for ε and 347.3 and 300 pm for σ, respectively. The adsorption isotherms and heats are given in Figure 2 and Table 2, respectively. The atomic positions in the systems lead to similar adsorption heat values and adsorption isotherms. For instance, the heats of adsorption change by 0.2 kJ.mol -1 at maximum. The other parameters have a more important weight. Switching from silicate to aluminophosphate oxygen Lennard-Jones parameters changes the isosteric heats of adsorption of N 2 by ca. 1.0 kJ.mol -1 . The point charge and ab initio electrostatic potentials change the heats of adsorption by ca. 0.6 kJ.mol -1 . The point charge electrostatic potentials show higher heats of adsorption than that obtained with ab initio electrostatic potentials.
For the eight possible combinations of systems, Lennard-Jones and electrostatic potential parameters reveal common trends and features when the N 2 Gibbs excess adsorption isotherm data are analyzed (Fig. 2). The N 2 excess adsorption increases with N 2 fugacity before to drop above ca. 40-50 bar. The values remain positive in the considered range (i.e., 10 -3 to 200 bar). As also pointed out by the isosteric heat of adsorption, the Lennard-Jones parameters are the dominant ones. This is clearly seen at 200 bar when all curves end up in two. The largest amount of N 2 is obtained for adsorption in materials that use aluminophosphate oxygen Lennard-Jones parameters. At low or moderate fugacity, the trend is also well marked: the amount of adsorbed N 2 is lower in systems defined with aluminophosphate oxygen Lennard-Jones parameters than with silicate ones. Although the electrostatic potential parameter is not the dominant in the case of the adsorption of N 2 , it plays a role and alters the adsorption isotherms (Fig. 2). The importance and influence of atomic charges on adsorption isotherms has been observed before in simulations [3,[27][28][29]. At high fugacity, there is practically no effect of the electrostatic potential as mentioned previously. At low and moderate N 2 fugacity, the amount of N 2 is always larger when a point charge electrostatic potential is used.

N 2 Adsorption Isotherms and Isosteric Heats in ZON-Type Materials
Now, we consider the adsorption of N 2 from fugacity of 10 -3 to 200 bar in materials 1, 2, 3 and 4. The Monte-Carlo simulations data are shown in Figure 3. The Gibbs excess sorption of N 2 increases in all solids before to reach a maximum at fugacity between 50 and 200 bar. The maxima are reached at 50 bar in 1 and ca. 75 bar in 2. When defects are present, the maxima are reached at higher fugacities: the maxima are reached at 100 and 200 bar in 4 and 3, respectively. The N 2 fugacity at which the Gibbs excess sorption maxima are reached follows the same ordering than that of the ZON porosity or micropore volume (Tab. 1). If we project the micropore volume into a hypothetical single cylinder contained in the solid unit cell, the diameters for 1, 2, 3 and 4 are 910, 760, 590 and 670 pm, respectively. To a large extent, the isosteric heats of adsorption of N 2 are constants in a given solid. It is only at N 2 fugacity close to that of the fluid condensation that they, first, increase slightly before a large drop. The small maxima in the adsorption heats do not exactly take place at the same fugacity than that of the Gibbs excess sorption (Fig. 3). The heats between 10 -3 to ca. 1 bar are approximately constant and equal to 16.8, 17.7, 18.4 and 18.8 kJ.mol -1 for 1, 2, 3 and 4, respectively. To compute the Gibbs excess sorption, we performed the Grand Canonical simulation of N 2 alone. This simulation gives us the isosteric heat of N 2 as function of fugacity: it is equal to 2.5 kJ.mol -1 up to ca. 10 bar and then increases up to 4.5 kJ.mol -1 at 200 bar. This does not provide support to better interpret or analyze the evolution of the isosteric heat of adsorption of N 2 in 1, 2, 3 and 4 as a function of the fugacity. The zero-pressure extrapolation heat of adsorption of N 2 in silicalite is experimentally reported to be 17.6 kJ.mol -1 [30]. We did not sample higher pressure than 200 bar but the behavior observed for 1 with N 2 fugacity beyond the adsorption heat maximum is likely to be equivalent for 2, 3 and 4.

CO 2 Adsorption Isotherms and Isosteric Heats in ZON-Type Materials
The Gibbs excess and isosteric heat of adsorption of CO 2 in the ZON materials obtained from the Monte-Carlo simulation are reported in Figure 4. The Gibbs excess sorption isotherms of CO 2 in 1, 2, 3 and 4 increase before they drop between 25 and 50 bar when CO 2 becomes liquid. The experimental partial pressure of CO 2 is 64.3 bar at 298 K [31]. The maximum value that is reached in 1 is larger than those in 2, 3 and 4, which are all three rather close. After condensation of CO 2 took place in the pores above 50 bar, the adsorption curves show equivalent slopes although they are shifted with respect to each other. At 200 bar, the adsorbed CO 2 amount is -0.05 g.cm -3 in 1 and 2, 0.04 g.cm -3 in 3 and -0.01 g.cm -3 in 4. We did not investigate adsorption above 200 bar but the curve slopes suggest that Gibbs excess sorption becomes negative in all solids.
The isosteric heats of adsorption of CO 2 in the different solids are again practically constant when the CO 2 fugacity is below 1 bar. The heats are 25.4, 27.0, 30.9 and 28.6 kJ.mol -1 , in 1, 2, 3 and 4, respectively, while the experimental heat of adsorption in silicalite is 27.2 kJ.mol -1 [30]. Above 1 bar, the heats of adsorption in 1 and 4 increase up to the value in 3, which remains at similar value than that at low fugacity. The heat in 2 increases by 0.5 kJ.mol -1 between 5 and 20 bar and, then, drops down to 24.7 kJ.mol -1 at 200 bar. The isosteric heat of pure CO 2 in the absence of solid is practically constant up to 30 bar and equal to 2.5 kJ.mol -1 and then it increases up to 27 kJ.mol -1 at 200 bar. This does not help to interpret or analyze the evolution of the isosteric heat of adsorption of CO 2 in 1, 2, 3 and 4 as a function of the fugacity.
The isosteric heats of adsorption of CO 2 and N 2 and especially their variations as a function of the fluid fugacity are qualitatively distinct. The heat of adsorption is practically unchanged for N 2 in the range 10 -3 -200 bar in 2, 3 and 4 and shows a drop above 50 bar in 1. In the case of CO 2 , only the adsorption heat in 3 remains practically unchanged. The heat of adsorption in 1 and 4 increases and it is in 2 that it drops at high fugacity.

H 2 O Adsorption Isotherms in ZON-Type Materials
AlPO 4 and SiO 2 microporous crystals are hydrophobic materials [32], the latter being less hydrophobic than the former. However, less is known when defects are present, although ion-exchanged zeolites and aluminophosphate crystals show increasing hydrophilicity with the amount of ions in the structure [32]. We analyze the water adsorption isotherms in our ZON models by means of Grand Canonical Monte-Carlo simulations. The results are shown in Figure 5: we report the total amount of adsorbed water per volume of solid in the inset graph, while the data are expressed per micropore volume of solid in the main graph.
The ZON micropore saturation with water occurs at higher fugacity than in the case of water alone, which takes place at 27 mbar for our water model. In our solids, the condensation of water is obtained between 0. Excess adsorption isotherms and isosteric heats of adsorption (inset) of N 2 in ZON-type materials. 1 and 10 bar in 2, 3 and 4. This is well above the water saturation pressure at 298 K. Hence, "wet" gases should not lead to the water filling of the solid micropores here. However, 4 shows a different water isotherm than the other solids. Even at extremely low water fugacity (1 mbar), a water molecule binds to the solid: inspection of the conformations generated by Monte-Carlo simulations reveals that this water molecule is always found in the vicinity of Na + . Unless the driest conditions are met, it seems experimentally impossible to avoid presence of this water molecule and the stoichiometric formula of 4 is better considered as [H 2 O•NaAl 32 SiP 31 O 128 ]. The amount of water in 1, 2 and 3 remains negligible below ca. 0.1 bar.

Adsorption Isotherms of "Wet" CO 2 -N 2 Mixture in ZON-Type Materials
Now, we focus on the results obtained on the Monte-Carlo simulations of the adsorption of CO 2 -N 2 mixtures (Fig. 6).
We analyzed the water adsorption and observed that water is present in 4 and can be considered absent in 1, 2 and 3.
For the adsorption of CO 2 -N 2 in 4, we set a "residual" H 2 O fugacity of 1 mbar and obtained in all simulation a water molecule present near Na + . We determine the CO 2 /N 2 molecular ratio as a function of fugacity in an empty cell by Grand Canonical Monte-Carlo simulations. The values are reported in the inset graph in Figure 6. Between 1 and ca. 30 bar, the CO 2 /N 2 ratio remains constant and equals to 1: CO 2 and N 2 behave like ideal gases. Between 50 and 80 bar, CO 2 condensates and the CO 2 /N 2 ratio becomes larger than 1. The values in the main graph in Figure 6 are normalized with the ratio values of the inset graph. The fugacity values in the x-axis are the sum of the CO 2 and N 2 fugacities, which are equal. There are more than ten times more CO 2 molecules adsorbed in the ZON micropores than N 2 molecules. This can be explained by the larger isosteric heat of adsorption of CO 2 than that of N 2 in the solids (Fig. 3, 4). In all solids, the CO 2 /N 2 ratio drops gradually above 10 bar and the ratios are around 3 in 2, 3, 4 and 5 in 1 at 100 bar. For all systems except 1, the ratios have an increasingly negative slope as a function of the fugacity. The case of 1 is different: the CO 2 /N 2 ratio actually increases between 1 and ca. 30 bar and drops at higher fugacity. The molecular ratios and their evolution as a function of the fugacity are practically equivalent in 2 and 4. The material that shows the highest affinity for CO 2 at 1 bar is 3, with a CO 2 /N 2 ratio of 15, while it is 13 in 1 and 12 in 2 and 4. Above 10 bar and because of the positive slope of CO 2 /N 2 ratio in 1, CO 2 shows the highest affinity in 1. The highest CO 2 affinity in 1 is reached at 30 bar and then it lowers when CO 2 starts to condensate in the micropores. The CO 2 and N 2 adsorption competition in the ZON micropores is difficult to estimate a priori from the adsorption of the pure fluids. The adsorption isotherm simulations of the binary (ternary with H 2 O) mixture appear mandatory to get an estimate of the CO 2 versus N 2 adsorption affinity in the solids. In the current work, we did not analyze time-dependent properties, which are also required to get a complete estimation of the separation properties of fluids by solids [33,34]. The Monte-Carlo formalism is not the most appropriate to simulate these and dynamic molecular simulations would be required to complement the evaluation of the separation properties of the above materials. However, due to the large similarities between the solid structures, it is likely that the diffusion constants in the different materials should not be very different.

CONCLUSION
We have shown that an intimate coupling of periodic DFT and MM Monte-Carlo methods is possible and permits to describe complex systems. Integration of periodic DFT and molecular mechanic methods in our chemical simulation interface MedeA permits to prepare and submit large number of calculations and to easily analyze them. The periodic DFT calculations give rigorous and unambiguous parameters that can be used directly in the molecular simulations. The bottleneck of the availability of forcefield parameters for the simulations is strongly reduced. The application of this approach to the separation of CO 2 /N 2 mixtures in presence of water by solids illustrates its sensitivity and elegance and its high chemical engineering interest. The simulations reveal that ZON with chemical composition Al 15 Si 2 P 15 O 64 is the material that shows the largest affinity for CO 2 adsorption below 10 bar. At higher fugacity, ZON with chemical composition SiO 2 becomes more suitable material to adsorb favorably CO 2 .