Molecular simulation of gas adsorption and diﬀusion in a breathing MOF using a rigid force field †

a Simulation of gas adsorption in flexible porous materials is still limited by the slow progress in the development of flexible force fields. Moreover, the high computational cost of such flexible force fields may be a drawback even when they are fully developed. In this work, molecular simulations of gas adsorption and diﬀusion of carbon dioxide and methane in NH 2 -MIL-53(Al) are carried out using a linear combination of two crystallographic structures with rigid force fields. Once the interactions of carbon dioxide molecules and the bridging hydroxyls groups of the framework are optimized, an excellent match is found for simulations and experimental data for the adsorption of methane and carbon dioxide, including the stepwise uptake due to the breathing eﬀect. In addition, diﬀusivities of pure components are calculated. The pore expansion by the breathing eﬀect influences the self-diﬀusion mechanism and much higher diﬀusivities are observed at relatively high adsorbate loadings. This work demonstrates that using a rigid force field combined with a minimum number of experiments, reproduces adsorption and simulates diﬀusion of carbon dioxide and methane in the flexible metal–organic framework NH 2 -MIL-53(Al).


A Introduction
Metal-organic frameworks (MOFs) consist of metal or metal clusters (''connectors'') linked by organic ligands (''linkers'') [1][2][3][4][5][6] forming nanoporous crystalline solid materials. The ligands act as spacers, creating an open porous structure with very high pore volume and surface area. Due to their unusual variety in terms of chemical composition, accessibility and pore dimensions as well as to their low densities (0.2-1 g cm À3 ) and high surface areas (500-4500 m 2 g À1 ), MOFs are considered promising candidates in applications involving adsorption and separation of strategic gases such as carbon dioxide, methane, and hydrogen. 3,4,[6][7][8][9][10] The first report on gas adsorption using MOFs was published in 1997 11 and since then, MOFs have become an active field of research, resulting in numerous publications. [12][13][14][15][16] In case of carbon dioxide adsorption, several approaches have been followed to improve selectivity: adsorbate-surface interactions, flexibility/gate opening mechanisms, and cooperative effects 14 (i.e. flexibility and specific interactions) have been reported as very effective ways of improving selectivity. In view of the high selectivities, relatively low enthalpies of adsorption, and fast kinetics, MOFs are among the most promising solids for carbon dioxide capture and separation.
Several authors have applied molecular simulations to study the adsorption and separation of carbon dioxide in metalorganic frameworks, i.e. Wang et al., 17 Babarao et al., 18 Keskin and Sholl, 19 Walton et al., 20 Yang and Zhong, [21][22][23] Calero et al., 24 Chen and Jiang, 25 Liu et al., 26 and Wells and Chaffee. 27 Also combined experiment and modeling techniques are used in work by Pera-Titus et al., 28 Lescouet et al., 29 Serra-Crespo et al., 30 Chen et al., 25 Stavitski et al., 31 Boutin et al., 32 and Montoro et al. 33 on functionalized MOFs. Although the origin of the improved separation performance in functionalized MOFs is not yet fully understood, it is clear from the above mentioned works that a delicate interplay of weak dispersion forces with framework polarity plays an important role.
A special class of MOFs is formed by those whose pore dimensions may change without breaking chemical bonds within the framework. This results in special properties like the breathing effect 8,34 and the gate phenomenon 35,36 where pores contract or open during molecule adsorption. An example of a breathing type material is the MIL-53 series. MIL-53 is built from MO 6 octahedra (where M can be three-valent ions Fe 3+ , Cr 3+ , Al 3+ , Ga 3+ , In 3+ or Sc 3+ ) formed from trans bridging OH ions and the oxygens of coordinate, bridging 1,4-benzenedicarboxylate linkers. 34 In this way a crystalline material with 1-D diamond shaped pores is formed. Amino-MIL-53(Al) 37 is a material with the topology of MIL-53. During the synthesis of amino-MIL-53(Al), 2-amino terephthalic acid is used as the linker molecule instead of terephthalic acid. The isoreticular material obtained has the following formula: Al(OH)[O 2 C-C 6 H 3 NH 2 -CO 2 ].
During the last few years, we have extensively worked on NH 2 -MIL-53(Al), demonstrating the high potential of this functionalized MOF in carbon dioxide separation and its thrilling optical properties. 16,31,32,[38][39][40][41] Both, the high carbon dioxide selectivity and the optical switch behavior of this structure are due to its unique breathing behavior. In the presence of carbon dioxide, two different pore configurations were found. The narrow-pore (np) has a smaller cell and pore volume and it is the dominant phase at lower pressures. Once certain pressure is reached, there is a pore expansion to the large-pore configuration (lp). 16 The chemical composition remains the same, and these configurations only differ in their pore dimensions.
The practical implementation of transition between phases has proved challenging from a simulation point of view. The development of specific force fields for every different structure and the intensive computational efforts needed to run simulations using flexible force fields certainly play against the goal of simulations taking over experimental efforts for predictive purposes. [42][43][44][45] In this work, we demonstrate that by using a rigid force field and with a minimum number of experiments, it is possible to reproduce adsorption and to simulate diffusion of carbon dioxide and methane in the flexible metal-organic framework NH 2 -MIL-53(Al).
The remainder of this paper is structured as follows. In Section B, the details of the simulation methodology and experiments are presented, including descriptions of the simulation techniques, force fields, and models used for the molecules and adsorbents. In Section C, the computed adsorption isotherms for carbon dioxide, Henry coefficients, and enthalpies of adsorption at zero coverage of both carbon dioxide and methane, the adsorption of the binary CO 2 :CH 4 mixture, and the diffusion of carbon dioxide molecules inside the pores are analyzed. Finally we give some concluding remarks.

B Simulation and experimental section
Atomic charges are used for all the atoms in the structure. In order to calculate the atomic charges for the narrow and largepore structures of NH 2 -MIL-53(Al), we used the Mulliken charge partitioning scheme 46 from all electron GGA DFT calculations employing the PW91 functional 47 and a double numerical plus polarisation (dnp) basis set. 48 The calculations were carried out with the DMol3 code. 49 We did not use the cluster approach to calculate the charges, because breaking the bonds can lead to spurious effects. We took into account the whole unit cell in the process of calculation of the charges. Ramsahye et al. 50 showed that the charges obtained with this method are less dependent on the geometry than charges obtained with electrostatic potential fitting methods of cluster models of the periodic structures, and they employed these charges to model successfully the adsorption of various molecules in MIL-47 and MIL-53. [50][51][52] Guest molecules are considered as small rigid molecules. Carbon dioxide is modeled as a linear molecule with a C-O bond-length of 1.149 Å, and partial charges are distributed along each molecule to reproduce the experimental quadrupole moment. 53,54 Methane molecules are described with a united atom model, in which each molecule is treated as a single interaction centre. 55 The energy of each gas atomic interactions is described by Lennard-Jones and Coulombic potentials. Lorentz-Berthelot mixing rules were used to calculate mixed Lennard-Jones parameters.
Monte Carlo (MC) simulations were performed to study the adsorption of carbon dioxide and methane in the Grand Canonical ensemble. To compute the Henry coefficients and enthalpies of adsorption at zero coverage of the studied molecules simulations were performed in the NVT ensemble using the narrow-pore structure. Simulations were performed using the in-house RASPA code. 56 The simulations were performed in cycles, and in each cycle, a MC move was chosen at random with a fixed probability: translation, regrowth, rotation, and insertion/deletion of molecules (in the Grand Canonical ensemble). We used at least 10 7 cycles and charge interactions were computed using Ewald sums with a relative precision of 10 À6 . Absolute adsorption was converted to excess adsorption 57 for comparison with experimental values.
Molecular dynamics (MD) simulations were carried out to compute the self-diffusivities of the carbon dioxide molecules inside the pores. The self-diffusivities were obtained by calculating the slope of the mean-square displacement at long times. 58 In these simulations the velocity-Verlet algorithm was used to integrate Newton's law of motion. MD simulations of 10 7 cycles in the NVT ensemble (Nóse-Hoover chain thermostat) were performed using a time step of 0.5 fs. 10 4 equilibration cycles were used and we took the initial positions of the molecules from previous MC simulations. The Lennard-Jones potential is cut and shifted with the cutoff distance set to 12 Å. NH 2 -MIL-53(Al) structures are modeled as rigid structures with Lennard-Jones parameters taken from the universal DREIDING force field. 59 Simulations were performed at 257, 263, 273, and 298 K using 32 unit cells (a = 19.7939 Å, b = 7.8496 Å, c = 6.5934 Å, a = g = 901, b = 73.56731) for the narrow-pore structure and 16 unit cells for the large-pore structure (a = 17.4426 Å, b = 12.0429 Å, c = 6.8447 Å, a = g = 901, b = 89.991). The crystal structures previously reported by our group 16 were used for the simulations. The pore volumes of the structures during the simulations are 0.07 cc g À1 and 0.52 cc g À1 , for the np and lp forms respectively.
Adsorption isotherms of carbon dioxide (purity of 99.995%) and methane (purity of 99.95%) were determined using the volumetric technique with an apparatus from BEL Japan (Belsorp HP). Around 0.5 g of NH 2 -MIL-53(Al) was placed in the sample container. Before every measurement, the adsorbent was pre-treated by increasing the temperature to 473 K at a rate of 10 K min À1 under vacuum and maintaining the temperature for two hours. The measurements for carbon dioxide and methane adsorption were carried out at 257, 263, 273, and 298 K.

C Results
The experimental adsorption isotherm for carbon dioxide in NH 2 -MIL-53(Al) at 257 K was compared with the adsorption isotherms obtained from simulation using both narrow-pore and large-pore crystallographic structures. The experimental isotherm shows a pronounced step, attributed to the transition in the framework structure due to the breathing mechanism. Carbon dioxide adsorption isotherms display the typical s-shape of breathing materials like regular MIL-53. The first plateau is reached at 3 mol kg À1 and pressures below 700 kPa. At higher pressures a drastic increase in the amount adsorbed occurred, results in a saturation loading of 10 mol kg À1 . The same behavior is found over a wide range of temperatures. The experimental saturation loading (B10 mol kg À1 ) is reproduced using simulations for the large pore structure. However, experimental adsorption in the low pressure range is underestimated. In order to address this issue, we increased the cell volume of the np configuration 16 by a 0.65% (increase of the cell volume of 0.65%), resulting in a slight increase in loading in the simulated adsorption isotherm.
In order to obtain a better match between experimental and simulated data in the np configuration, the framework-adsorbate interactions were optimized using a universal force field. Although universal force fields are usually expected to provide realistic predictions of molecular structures, when dealing with complex adsorbents like MOFs, these generic force fields usually deliver poor results. Previous studies pointed at a strong interaction between carbon dioxide molecules and the framework hydroxyls. 60 Using experimental data, we could quantify the extent of this interaction by direct fitting of the data on the isotherm ( Fig. 1(a)). C CO 2 -H OH Lennard-Jones parameters were obtained by calibrating the force field through explicitly fitting one point (10 kPa) of the experimental isotherm at 257 K. After calculating the appropriate Lennard-Jones interaction at this single temperature, the full isotherm could be reproduced and predicted over a wide range of temperatures (see Fig. 1 (b) 263 K, (c) 273 K, and (d) 298 K) using the np structure, resulting in an excellent agreement with experimental data (full green symbols).
Simulated zero-coverage adsorption enthalpies and Henry coefficients for carbon dioxide and methane are compared with previous experimental data 39 in Table 1. Simulated adsorption values are very close to the experimental data. Results for carbon dioxide confirm a strong interaction between the molecule and the structure. In case of methane, a much larger difference between experimental and simulated K H is evident. We attribute this difference to the large error in the calculation of the experimental value due to the very low affinity of the framework for methane.

View Article Online
The main goal of this work is to capture the nplp transition in NH 2 -MIL-53(Al) at a low computational cost. In Fig. 2, experimental adsorption isotherms of carbon dioxide at different temperatures are shown as a function of the relative pressures (p/p 0 ). The adsorption isotherms present a similar behavior. In the studied temperature range, the nplp phase transition occurs at the same relative pressure of about p/p 0 = 0.3, where p 0 is the pressure of the vapour in equilibrium with its non-vapour phases. This fact clearly infers that the phase transition is related to the chemical potential of the adsorbate. Based on our previous crystallographic results, 16 a first estimation could be made of the relative population of the np and lp phases during framework expansion. Based on that estimation of the relative population of each phase, a linear combination of the simulated isotherms using the rigid structures as function of p/p 0 was constructed. An excellent match between the isotherm obtained by the linear combination of the rigid simulated isotherms and the experimental data was obtained ( Fig. 3(a), blue asterisks). By extrapolating the obtained results to higher temperatures, the transition between the two crystallographic phases of NH 2 -MIL-53(Al) can be well reproduced at different temperatures from only one experimental isotherm and a linear combination of the percentages of each coexisting structure at a given p/p 0 . Following this approach, both the experimental and the computational costs are minimized, and flexible force fields are not needed to capture the breathing effects due to the interaction with guest molecules. Fig. 3 shows the experimental adsorption isotherms (black line) at (b) 263 K, (c) 273 K, and (d) 298 K compared with the results obtained via linear combination of the two coexisting structures (blue asterisks). For the whole temperature range studied, agreement between simulation and experimental data is excellent.
Adsorption isotherms of individual components and equimolar carbon dioxide and methane mixtures as a function of the partial pressures at 298 K computed using MC simulations are shown in Fig. 4. The adsorption of (top) carbon dioxide and (bottom) methane is presented for the single component (full symbols) and as a part of the binary mixture (empty symbols) in the narrow pore form (circles and diamonds, respectively), in the large pore form (squares and triangles, respectively) and in the linear combination of the two coexisting structures (asterisks). Carbon dioxide starts adsorbing at lower pressures than methane as single component and has a 33% higher  saturation loading. The adsorption of carbon dioxide is almost independent of the presence of methane.
In contrast, when methane is part of the binary equimolar mixture, its equilibrium loading is much smaller than when it adsorbs as single component. The strong decrease in methane adsorption in the presence of carbon dioxide is due to two different reasons: (i) on one hand, carbon dioxide molecules are already confined in the pores before methane can be adsorbed due to the constrained pore dimensions of the np form and (ii) the different interaction strengths between the two adsorbates and the hydroxyl groups of the metal-organic framework. Due to the significant quadrupole moment of carbon dioxide, the interaction with the bridging hydroxyls is much higher than for the non-polar methane molecule. As a result, NH 2 -MIL-53(Al) is a highly selective adsorbent.
MD simulations are a powerful tool to obtain information concerning gas diffusion for single components in MOFs. Up to now, only a few computational studies based on MD simulations have been reported so far on the diffusion of carbon dioxide in MOFs. [61][62][63][64][65] As we have recently shown, 38 the interaction between amino groups of the linker and the hydroxyl groups of the structure prevents rotation of the linker, in contrast to the unfunctionalized MIL-53 structure. Therefore, diffusion of molecules from one channel to another is highly unlikely, making directional diffusion along the channels the most likely diffusion mechanism in NH 2 -MIL-53(Al).
In Fig. 5, calculated self-diffusivities as a function of the number of adsorbed molecules per unit cell at 273 K are presented. The very slow diffusivity for less than 3 molecules per unit cell, the maximum loading for the close configuration of NH 2 -MIL-53(Al), is related to the fact that molecules are not able to pass each other due to their orientation into the one dimensional channels. Diffusivity is faster once the structure is open, even when a larger number of molecules are adsorbed per unit cell. The latter is due to the larger available space in the pore and the weaker interaction with the pore walls. Fig. 6 shows two snapshots of carbon dioxide in the np (left) and lp (right) structures of NH 2 -MIL-53(Al) taken directly from our simulations at 273 K. Although we are not able to provide a comparison with experimental data on NH 2 -MIL-53(Al), on the basis of the excellent agreement found between experiments and simulations for equilibrium adsorption data, we are confident that the diffusivities here reported are qualitatively correct, as it has been shown previously for zeolites using similar molecular simulation methods. 66,67

D Conclusions
Adsorption of carbon dioxide and methane in the flexible metal-organic framework NH 2 -MIL-53(Al) has been simulated, both for pure component and equimolar binary mixtures, using rigid force fields for the framework. The use of flexible force fields was avoided following a two-steps process, with which we   obtained a very good agreement between the available experimental data and simulations at different temperatures. The first step is the optimization of the interactions of carbon dioxide molecules and the bridging hydroxyls groups of the framework followed by the use of a linear combination of the np and lp forms based on experimental results. The step found in the adsorption isotherms of carbon dioxide due to the breathing effect is well reproduced by our simulations, and the high selectivity towards carbon dioxide in the presence of methane is explained on the basis of size exclusion and affinity. Computed enthalpies of adsorption and Henry coefficients for both gases are also in agreement with experimental data. Computed selfdiffusivities of carbon dioxide in NH 2 -MIL-53(Al) show two different regions. The pore expansion by the breathing effect influences the self-diffusion mechanism and much higher diffusivities are observed at relatively high adsorbate loadings.
Summarizing, in this work we have shown that adsorption and diffusion of gases can be simulated in flexible porous materials without employing flexible force fields and with very limited experimental efforts.