Adsorption of Light Gases in Covalent Organic Frameworks: Comparison of Classical Density Functional Theory and Grand Canonical Monte Carlo Simulations

A classical density functional theory (cDFT) based on the PC-SAFT equation of state is proposed for the calculation of adsorption equilibria of pure substances and their mixtures in covalent organic frameworks (COFs). Adsorption isotherms of methane, ethane, n-butane and nitrogen in the COFs TpPa-1 and 2,3-DhaTph are calculated and compared to results from grand canonical Monte Carlo (GCMC) simulations. Mixture adsorption is investigated for the methane/ethane and methane/n-butane binary systems. Excellent agreement between PC-SAFT DFT and GCMC is obtained for all adsorption isotherms up to pressures of 50 bar. The cDFT formalism accurately predicts the selective accumulation of longer hydrocarbons for binary mixtures in the considered COFs. This application shows substantial predictive power of PC-SAFT DFT solved in three-dimensional geometries and the results suggest the method can in the future also be applied for efficient optimization of force field parameters or of structural properties of the porous material based on an analytical theory as opposed to a stochastic simulation.


Introduction
case of COFs, containing significantly larger pore sizes compared to zeolites, the impact of the stacking motifs of adjacent layers in the structural model used in the simulations is expected to influence the simulation outcome at least to the same extent as residual discrepancies in the force fields used. 21,[26][27][28] To increase versatility of computational methods, a more efficient alternative to molecular simulation could be classical density functional theory (cDFT) 30,31 which is also rooted in the framework of statistical mechanics but relies on an inhomogeneous density profile compared to explicit atomistic molecular simulation. One of the most common application of cDFT is adsorption in homogeneous slit pores with two opposing planar walls. The solid is thereby modeled by an external field which commonly takes the form of a Lennard-Jones 9-3 potential or a Steele potential. cDFT accurately predicts the adsorption behaviour when compared to GCMC simulations including layering transitions 32 and capillary condensation. 33 The adsorption behaviour of real unordered porous materials, however, is often not well represented by the homogeneous slit pore model with one given pore size. This is because of the ambiguous pore structures with often unknown porosity, chemical composition and pore size distributions. Therefore, cDFT models include heterogeneities, 34 both in pore size distribution and surface roughness/chemical heterogeneity, to compare accurately to adsorption experiments. Ordered porous media, in turn, are characterized by their regular pore structure and, thus, provide a consistency test between cDFT and molecular simulations beyond one-dimensional homogeneous slit pores.
A key ingredient of cDFT is the Helmholtz energy functional used to describe the fluid-fluid interactions. Whereas the hard-sphere repulsion is often represented by a functional based on Fundamental Measure Theory, [35][36][37] dispersive attractions are either treated by a simple mean-field theory which ignores density correlations of the fluid or non-local weighted density approximations in combination with an underlying equation of state. Comparative computational studies of the adsorption in ordered porous frameworks between cDFT and GCMC simulations were performed by different groups; each utilizing different Helmholtz energy functionals. Guo and coworkers 38 compared adsorption isotherms of noble gases in metal organic frameworks using a mean-field approach. Fu and Wu 39 assessed the performance of different dispersive Helmholtz energy functionals from mean-field theory to weighted density approximations with an empirical equation of state. They compared calculated adsorption isotherms of methane in metal organic frameworks to molecular simulations.
This study uses a functional based on the Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) equation state, 40,41 which also utilizes a weighted density approximation. 42 This functional was already successfully applied to adsorption in one-dimensional slit pores 33 and the calculation of surface tension and Tolman lengths. 43 In this study we assess the PC-SAFT DFT model for predicting adsorption in ordered threedimensional COF frameworks. We consider the adsorption of light gases in two typical COFs and we compare results from GCMC and cDFT. The results are discussed in light of methodological differences of the two approaches.

Classical Density Functional Theory
In this section, we summarize the fundamental equations of classical density functional theory and the application to adsorption in COFs. Density functional theory is formulated in the grand canonical ensemble at constant chemical potentials = , = 1, … , of all species, volume , and temperature . The grand canonical potential was shown to be a unique functional of the inhomogeneous density profile ( ) = ( ), = 1, … , and can be expressed as where [ ( )] is the intrinsic Helmholtz energy functional capturing the fluid-fluid interactions and ext ( ) is the external potential due to solid-fluid interactions acting on species . For adsorption in microporous materials it is instructive to think of the system as being connected to a large bulk reservoir with the same temperature and chemical potentials , so that a pressure of a communicating bulk fluid ( , ) can be calculated. The equilibrium density distribution 0 ( ) minimizes the grand canonical functional and its value is then equal to the grand canonical potential Ω ( , , ), so that The equilibrium density profile is obtained by solving the Euler-Lagrange equation using a damped Picard iteration in combination with an Anderson mixing scheme to accelerate the convergence rate. 44 For a compact notation, we henceforth omit the subscript 0 in the equilibrium density profile; we use ( ) for the vector of density profiles of all components in the system. The with repulsive hard-sphere interactions 36,37,45 (hs), hard-chain formation 46,47 (hc), and van der Waals leading to ( ) = ( ) for homosegmented chains.
The external potential represents the van der Waals interactions exerted by the COF atoms onto a fluid (segment). The external potential is calculated by considering the interactions of a PC-SAFT molecule with all individual solid atoms of the framework, leading to where is the number of solid atom interaction sites of the considered framework and is the where and are the Lennard-Jones interaction parameters of atom interaction site taken from the DREIDING force field 50 representing the COF structure.
In this work, the vector containing the number of adsorbed molecules ads = ads , = 1, … , of a component mixture is calculated with for determining the equilibrium densities ( ) and subsequently obtain the adsorbed amount using eq. (8).

Grand Canonical Monte Carlo Simulation
All GCMC calculations were performed using the molecular simulation software RASPA. 51 Intramolecular fluid and intermolecular fluid-fluid interactions were described with the TraPPE force field. 52,53 The CH x groups in methane, ethane and -butane were considered as single, chargeless interaction centers (united atoms) with effective Lennard-Jones potentials. Parameters for unlike interaction sites were determined using Lorentz-Berthelot combining rules. TraPPE approximates the quadrupolar nature of nitrogen by placing negative partial atomic charges at the position of the nitrogen atoms and a neutralizing positive partial charge at the center of mass. The COF framework was considered to be rigid such that only Lennard-Jones parameters and partial atomic charges needed to be assigned to the different atomic species. The Lennard-Jones parameters were taken from the DREIDING force field. 50 Partial atomic charges of the COF structures were calculated using the extended charge equilibration (EQeq) 54 41 Although the density beyond the cut-off radius is not uniform, we applied analytic corrections to the long-range Lennard-Jones tail, in order to reduce the sensitivity of the results with respect to the cut-off radius. 61 The real part of the electrostatic interactions was evaluated up to a cut-off radius of 12.0 Å. Long-range electrostatic interactions were calculated by Ewald summation 62,63 with a relative precision of 10 −6 . To carry out simulations at constant chemical potential, the PC-SAFT EoS was used to pre-compute a fugacity coefficient at the given temperature and pressure that was then passed to the MC code. The number of MC cycles was 25 × 10 4 , both, for equilibration and for the production phase. One cycle consists of max(20, ) MC moves (with as the sum of adsorbate molecules in the system), i.e. translation, insertion or deletion and, in case of molecules represented by more than one site, rotation moves.
In simulations of binary mixtures identity swap moves were carried out additionally. All moves were performed with equal probability.

COF Structures
The two COFs considered in the present work are the ketoenamine-linked COF TpPa-1 64 and the imine-linked COF 2,3-DhaTph 65,66 having pore sizes of approximately 1.8 and 2.0 nm, respectively, see Figure 1.
As the purpose of this study is the comparison between two computational approaches some- For the COF 2,3-DhaTph initial coordinates in a perfectly eclipsed arrangement were taken from the CoRe COF database. 12,13 However, the layer-layer distance in that structure of 6.7 Å is much larger than the experimentally reported value of 4.0 Å because the benzene rings were rotated by 90 • . Moreover, the lattice was not tetragonal as in the experimentally derived X-ray structure, 65

Ideal Adsorbed Solution Theory
Adsorption isotherms of mixtures can be estimated from the pure component isotherms using the ideal adsorbed solution theory (IAST). 68 In the present work the IAST equations were solved with the pyIAST package. 69 To account for non-ideal behavior of the gas phase at elevated pressure fugacities instead of pressures were employed in the IAST equations. 70,71

Results and Discussion
Before comparing adsorption isotherms predicted by cDFT and GCMC we first investigate vaporliquid equilibria to assess whether the two approaches show deviations that may impact their comparability. Note that the segment size parameter and the dispersive energy parameter used in PC-SAFT are different from the force-field parameters used in the MC simulations, even for methane. Both were independently adjusted to experimental data of pure compounds. Results from both approaches are comparable, however, because pure component parameters were adjusted to experimental data for phase equilibria.

Pure Component Adsorption
The pure component adsorption isotherms are presented by plotting the average absolute amount adsorbed ads per mass of COF-framework as function of the pressure in the external reservoir. The statistical uncertainties in the GCMC results are in almost all cases smaller than the symbol size. Figure 3 shows adsorption isotherms of nitrogen, methane and ethane in COF 2,3-DhaTph at 298 K with varying pressure up to 50 bar. For methane excellent agreement between cDFT and GCMC is obtained over the entire pressure range. For ethane the agreement between the two approaches is very good up to pressures of 2 bar. At higher pressures cDFT slightly underestimates the amount adsorbed. For nitrogen the agreement between GCMC and cDFT is remarkable given that the force field contains three collinear partial atomic charges to model the quadrupolar nature of the molecule while the PC-SAFT model entering the cDFT calculations describes nitrogen as non-polar, so that the van der Waals parameters effectively capture the (mild) quadrupole moment of nitrogen. Of course, the influence of the partial atomic charges on the adsorption behavior strongly depends on the magnitude of the charges of the adsorbent-framework, as reported for siliceous zeolites. 77 For here considered COF 2,3-DhaTph nitrogen adsorption isotherms with and without framework partial charges show only minor differences (see Figure S7 in the Supporting Information).  Chapman 78 and has also been applied with the PC-SAFT DFT model in previous work of our group. 79 For nitrogen the cDFT isotherm is slightly lower than the GCMC one. Again, we tested the influence of the partial charges on the framework atoms with respect to nitrogen adsorption and found that the GCMC isotherm is in excellent agreement with the cDFT isotherm when evaluated in a framework exempt from partial charges (see Supporting Information). The TpPa-1 framework is somewhat more polar than 2,3-DhaTph, if we use the sum of squared partial charges as a measure for how polar a framework is. We thus regard the sum of 2 ∕ , where is the number atoms of species in the simulation cell and its volume. For all of the four atomic species (C, H, N, O) these values are higher for TpPa-1 compared to 2,3-DhaTph. Therefore, for frameworks with low to moderate charge densities, adsorption of quadrupolar fluids may be approximated by dispersion interactions alone. In summary, and in view of the fact that no parameter is adjusted for relating the two modelling approaches, we consider the overall agreement observed in Fig. 3 and 4 as good.   relative to methane we introduced a second -axis in Fig. 6, to better present the amount of adsorbed methane. For methane, cDFT and GCMC are in very good agreement as can be expected from the comparison of the pure methane isotherms discussed above. The -butane isotherms as predicted from cDFT underestimate the adsorbed amount as compared to the GCMC results, similar to the behavior for pure -butane. For methane IAST predicts a slightly higher amount adsorbed above pressures of 10 4 Pa.

Conclusion
A classical DFT approach relying on a Helmholtz energy functional based on the PC-SAFT equation of state was used to predict adsorption equilibria of pure components and binary mixtures in covalent organic frameworks. The results were compared to adsorption isotherms from GCMC simulations (using the TraPPE force field for the fluids). While the latter approach is rooted in a fully atomistic description (within a united atom approximation for alkanes), cDFT employs a coarser The remarkable agreement of the two approaches, even for -butane, shows that cDFT is a powerful alternative to GCMC for studying adsorption equilibria in porous materials. The advantage of this approach is its analytical nature allowing the calculation of derivatives and, therefore, optimization tasks with respect to force field parameters or structural properties of the porous materials. The second aspect is in particular relevance for COFs because the stacking motifs have a substantial impact on the materials properties, including the adsorption behavior. The current limitation of the cDFT approach is a lower coverage of the chemical space compared to molecular simulations, in particular with regard to polar molecules.

Author Notes
Christopher Kessler and Johannes Eller contributed equally to this work.