Molecular dynamics study of CO₂ absorption and desorption in zinc imidazolate frameworks

Abstract

We report a study of the absorption of CO₂ into a number of zinc imidazolate framework structures, and subsequent desorption, using the molecular dynamics simulation method with force fields partly developed by ourselves. The simulations primarily give results concerning the mechanism of CO₂ absorption under conditions likely to be found in industrial waste gas streams. In particular we compare the rate of uptake of CO₂ for different ZIFs. We also show that it is possible to observe desorption by reduced pressure and high temperature. These characteristics confirm that ZIFs might be suitable for CO₂ absorption within an industrial capture and sequestration process.

---

Design, System, Application

ZIFs have been proposed as materials for CO₂ gas capture. We have developed a new simulation model and strategy to study CO₂ absorption in such materials, including a new model for CO₂ that works in both gas and condensed phases. This method can be adopted to study all kinds of gas absorption in ZIFs and other metal–organic-frameworks. It can also be used to study gas selectivity of MOF membranes. CO₂ absorption and de-absorption in a number of ZIFs have been studied. The results help us to understand the mechanism of the interaction between CO₂ and the ZIF structure, and to predict the absorption storage for similar structures. We show that ZIFs have great potential for CO₂ capture. For future application, a) our work shows that ZIFs have great potential for CO₂ capture; b) our method can be more widely used for many future applications. It can speed up studies of gas absorption, gas selection and storage by screening. This helps choosing appropriate MOFs for particular purposes such as toxic gas absorption in power plant, and hydrogen storage and transport for clean power, etc.

1 Introduction

The problem of the rapid increase of carbon dioxide in the atmosphere is well known, and in order to reduce releases from industrial processes capture and storage solutions have been proposed. One aspect of active consideration is to use porous crystalline materials to absorb CO₂ directly from gas streams. Metal–organic framework (MOF) materials^1,2 – porous crystalline materials formed as network structures where metal cations are linked together by organic ligands – have attracted a lot of interest in this regard. Their porous structures give the possibility of high molecular absorption capacity,^3,4 coupled with a high stability at different temperatures and pressures.

In this paper we study the absorption and desorption of CO₂ in one class of MOFs, namely the materials with the generic name zeolitic imidazolate frameworks (ZIFs).⁵ These have the chemical formula Zn(im)₂, where im = C₃N₂H₃⁻ or a related ligand. These materials form as networks of zinc cations linked by the im ligands. The angle subtended by the metal–im–metal vectors in many ZIFs (ca. 145°) is very similar to the Si–O–Si angle in silica, SiO₂, meaning that the networks are able to form structures whose topologies are either direct analogues of, or related to, siliceous zeolites. Furthermore, the structures formed by these materials can be tailored by changing the metal cation, tuning the exact imidazolate linkage ligands,⁶ and growing from different solvents. This has led to at least 100 types of ZIFs being found.^6,7 Their various porous structures offer the prospect of being exploited for applications in areas such as catalysis, storing gases such as H₂, and capture of gases such as CO₂.^6,8–13

Our approach is to use the molecular dynamics simulation method using empirical force fields, focussing on the ZIFs identified in the first report of these materials.⁶ In the work reported here we started with a pristine slab of a ZIF in contact with CO₂ gas, and applied an appropriate temperature and pressure. The pressure had the effect of pushing the CO₂ molecules into the ZIF slab. In this paper we report our observations on the processes by which the CO₂ can enter the ZIF slabs from the gas phase, and how the molecules move within the ZIF. Here we will also report the structures formed once maximum capacity has been more-or-less reached. We have performed simulations at different pressures, and we compare the results for the different ZIF systems. We have also investigated the desorption process induced by reduced pressure and increased temperature, but because of technical issues this proved to be rather harder.

In the next section we describe the methods used in this work, including details of the force fields we have developed, and the way we have used the molecular dynamics method in this work. Our simulation results are reported in the subsequent sections. The first set of results concerns the absorption of CO₂ into slabs of a range of different ZIFs, looking at both the time scales for absorption and then the CO₂ density profiles through the slab. Results are presented for two gas pressures. In order to understand the absorption process the section after considers in more detail the interactions between imidazolate ligands and CO₂ molecules. In the penultimate section we discuss the process of de-absorption with reduced pressure and increased temperature, with some conclusions presented in the final section. The ESI† presents more results than can be accommodated in the main paper.

2 Force fields

2.1 Ab initio calculations for parameter fitting

Because of the system sizes and time scales involved in our simulations, this work is only made feasible by using empirical force fields to describe the interactions between atoms. Some aspects of our force fields have been discussed previously.¹⁴ Briefly reviewing the ZIF model from our earlier work,¹⁴ we note that the force field we have used is a combination of new components based on ab initio calculations and previous highly-successful empirical inter-atomic potentials for organic molecules.¹⁵

Ab initio DFT calculations for small clusters of zinc cations and ligands were performed using the quantum chemistry program NWChem^16,17 with the standard PBE functional,¹⁸ and using either aug-cc-pVDZ or aug-cc-pVTZ basis sets depending on system size and following convergence tests. We used the distributed multipole analysis (DMA) method¹⁹ as performed by the program CamCASP²⁰ to calculate atomic multipole moments from the electronic wave functions, from which we obtained values for atomic point charges using program MULFIT.^21,22

The new interatomic potentials for our force fields were obtained by fitting to the energies for a large suite of atomic configurations obtained from ab initio calculations, also performed using the program NWChem,^16,17 but this time using the Møller–Plesset perturbation theory MP2 approximation. As for the DFT calculations, we used aug-cc-pVDZ and aug-cc-pVTZ basis sets. Values for the parameters in the force field equations (see immediately below) were obtained using the fitting procedure with the code GULP.²³

2.2 Basic equations in our force field

We used a number of different types of force-field functions to represent different facets of the set of interactions. For interactions between two bonded atoms we used the Morse potential:

E(r) = ε(exp(2α(r − r₀)) − 2 exp(α(r − r₀))) (1)

where r is the distance between the two atoms, and ε, α and r₀ are parameters of the model with values to be obtained by fitting to the ab initio energy surfaces. Bond-bending potentials for N–Zn–N and Zn–N–C angles were modelled using interactions of the form

E(θ) = (k/2)(θ − θ₀)² (2)

where θ is the bond angle, θ₀ is a parameter that represents the equilibrium angle (which we assumed to be a standard angle based on the local coordination, such as cos⁻¹ (−1/3) for the tetrahedral N–Zn–N angle), and k is a parameter of the model. We also included an inversion force in order to maintain the planarity of the group of atoms centred on the N atom together with the neighbouring Zn and C atoms. This has the form

E(ϕ) = A(1 − cos ϕ) (3)

where ϕ represents the angle between the normals to any two planes defined by a group of 3 out of the 4 atoms, and A is a parameter. In all these three potentials we excluded direct Coulomb interactions between bonded atoms.

For non-bonded interactions, such as between the atoms of different imidazolate ligands for example, we used the standard Buckingham form to represent both the long-range dispersion interaction and the short-range repulsion:

E(r_ij) = −(C_ij/r_ij⁶) + B_ij exp(−r_ij/ρ_ij) (4)

where r is the distance between atoms, i and j label the type of atom, and C_ij, B_ij and ρ_ij are parameters.

2.3 Force field models for the zinc imidazolate frameworks

The model for the interactions between the zinc atoms and the imidazolate molecular anions has been described in detail in our previous publication.¹⁴ In this work we treat the imidazolate ion as a rigid object with the atoms acting as interaction centres. Charges on both the imidazolate anion and zinc cation were assigned on the basis of DMA calculations of clusters. Similarly cluster calculations were used to parameterise the force fields involving the zinc cation and imidazolate anion. These included a Morse potential to describe the Zn–N bond (eqn (1)), a bond bending potential to describe flexing of the tetrahedral N–Zn–N and the Zn–N–C bond angles (eqn (2)), and an inversion potential (eqn (3)) to align the Zn–N bond within the plane of the imidazolate ligand. The development of these potential functions has been described in detail in our previous paper,¹⁴ where the parameters are listed.

For interactions between the imidazolate ligands we used the Buckingham potentials (eqn (4)), with parameters taken from the work of Williams.²⁴

2.4 Development and validation of a new force field model for CO₂ molecules in condensed and fluid phases

2.4.1 The model for interactions between CO₂ molecules. Whilst there are several models for the interactions between CO₂ molecules that are in widespread use in the literature, we preferred to develop a new model and tune the parameters to reproduce not only the gas-phase data but also crystalline data. We treated the molecule as a rigid body with a C–O distance of 1.16 Å.²⁵ The charges were obtained by the same DMA method described above, now using aug-cc-pVTZ basis sets allowed by the small size. The obtained values for charges are −0.30403e for C and 0.60806e for O. Interactions between CO₂ molecules were described using the atom–atom Buckingham potential (eqn (4)).

Values for the parameters in the Buckingham potential were obtained by fitting to a mixture of data. We used a suite of ab initio energies calculated for a large number of configurations of pairs of molecules of orientations as shown in Fig. 1. These ab initio energies were calculated using the MP2 method with aug-cc-pVTZ basis sets. We also made use of the crystal structure and phonon dispersion curves of solid CO₂ in the fitting procedure, which was carried out using the code GULP.²³ We also used Grand Canonical Monte Carlo (GCMC) simulations of the fluid phase diagram to help tune and constrain parameter models.

Fig. 1 CO₂ clusters used in MP2 calculations in support of the development of the force fields. The right-side molecules were moved horizontally to the right side to get sets of cluster energies with different distances between the molecules.

By this route we were able to tune values for all the parameters. The final set of parameters for our new CO₂ model are given in Table 1. Part of this work was briefly described in an earlier paper on fluid CO₂,²⁶ but here parameter values have been slightly refined.

Table 1 Fitted values of the parameters in Buckingham interaction (eqn (4)) between CO₂ molecules

Atom pair	B (kJ mol^−1)	ρ (Å)	C (kJ mol^−1 Å^−6)
C–O	190 909	0.2637	1216.21
O–O	203 567	0.2659	2149.19
C–C	108 347	0.2778	0.0

---

A model for CO₂ with non-rigid molecules are also developed and tested. A Morse potential was used for the stretching of the C–O bond, and a potential of form E(ϕ) = K(1 + cos θ) was used to describe flexing of the O–C–O angle θ; parameter values were tuned to give best agreement with the oft-reported frequencies of the internal molecular vibrations. The highest frequency vibration of the CO₂ molecule, namely the asymmetric C–O stretching mode, has a frequency of 70.5 Hz, which is much higher than CO₂ intermolecular vibration frequency (4 THz) as shown in Fig. 2, and higher than the highest network flexing frequency of the ZIFs of 14 THz.¹⁴ If we scale the timestep by the ratio of frequencies, 14/70.5 ≃ 0.2, we would require computational resources five times greater than we used for the rigid molecule model; given the very large cost of computer time incurred in this work, such an overhead would need significant justification. In the event our tests on the fluid phase diagram showed that the rigid and flexible molecule models gave the same results. Hence only the rigid body model was used here.

Fig. 2 Phonon dispersion curves of CO₂. The diagram on the left shows the calculation of our model (continuous curves) and the diagram on the right shows the corresponding calculations using the Potoff and Siepmann model.²⁷ The filled black circles are experiment data from inelastic neutron scattering.²⁹

2.4.2 Crystal structure and phonon dispersion curves. The CO₂ model has been validated by comparing with experiment data in the crystalline phase as well as the more-usual fluid phases. We also compare with calculations using the commonly-used model of Potoff and Siepmann.²⁷ In this model, the CO₂ molecules are treated as rigid bodies. Partial charges for O and C atoms are 0.7e and −0.35e respectively. Intermolecular interactions were modelling using the Lennard-Jones potential rather than the Buckingham potential preferred here. The parameters of Potoff and Siepmann²⁷ are reproduced in Table 2.

Table 2 Parameter for Lennard-Jones interaction used in the Pottoff and Siepmann model²⁷

Atom pair	ε (K)	σ (Å)
C–C	27.0	2.80
O–O	79.0	3.05
C–O	46.18	2.92

---

The crystal structure was taken from the X-ray refinement of Simon and Peters.²⁸ CO₂ has Pc3 space group symmetry, with the lattice parameter a = 5.624 ± 0.002 Å. The C–O distance is 1.15 Å. Our structure optimisation performed using GULP gave a lattice parameter 5.624 Å in perfect agreement with the experimental result, whereas the optimised lattice parameter for the Potoff and Siepmann model is 4.998 Å with 11% difference from the experiment. Typically we expect a good empirical model to have an accuracy for the lattice parameters of up to 2% (perfect agreement with experimental data should not be taken too seriously), and so we feel that the Potoff and Siepmann model has too poor level of agreement with experimental data in the crystalline phase.

To test the CO₂ model of the solid phase further, the phonon dispersion curves have been calculated using GULP²³ and compared with the inelastic neutron scattering data of Dolling et al.²⁹ (phonon measurements were performed for wave vectors along the [1,1,1] direction). Calculated phonon dispersion curves for both the model developed here and from that of Potoff and Siepmann are compared with the experimental data in Fig. 2. The calculations were performed with the relaxed structure. The experimental results were obtained by inelastic neutron scattering. From the figure we can see that our calculated dispersion curves show a good agreement with the experimental data. The results for the two acoustic modes and the range of optical phonon frequencies at wave vectors [0,0,0] and of our model are close to the experimental data. On the other hand, the phonon frequencies calculated using the Potoff and Siepmann model are systematically too high.

2.4.3 Grand canonical Monte Carlo simulation of CO₂ fluid. Our model for crystalline CO₂ has shown good consistency with experiments as discussed in the above section, a further study of the vapour–liquid phase has been simulated and compared with experimental data obtained from NIST database.³⁰ The simulations have been carried out using the GCMC code Towhee,³¹ which was written for vapour–liquid calculations.

GCMC simulations were performed for a range of temperatures from 200 K to 300 K, using the constant pressure and temperature Gibbs ensemble. The temperature and pressure coexistence lines are shown in Fig. 3 together with the density at the co-existence line. The simulation results show a satisfying agreement with the experiments. Detailed results for the vapour–liquid coexistence point are given in Table 3 and compared with results for the Potoff and Siepmann model. The accuracy of our model is comparable with that of Potoff and Siepmann for the fluid phase diagram, and has the advantage of also handling the crystalline phase much better.

Fig. 3 CO₂ vapour–liquid coexistence line. The black dot and star is the critical point from experiment³⁰ and our simulation respectively.

Table 3 Comparison of experimental³⁰ and calculated values for the CO₂ vapour–liquid coexistence point using the GCMC method. The two potentials are from the present work and the model of Potoff and Siepmann (PS)²⁷

	Temperature (K)	Pressure (bar)	Vapour density (g ml^−1)
Experiment	304.1	73.8	0.4676
Our model	308.4	80.8	0.4670
PS model	306.2	77.7	0.4649

---

2.5 Model of inter-molecular forces between CO₂ and the ZIF crystal

The values for the B and C parameters in the atom–atom Buckingham potentials (eqn (4)) used to describe the interaction between the CO₂ molecules and the Zn cations and ZIF ligands were obtained using again the procedure of fitting to a set of ab initio energies for a suite of configuration, as we described for the CO₂ molecules. The total energies of the system with different distances were calculated using the MP2 method with aug-cc-pVDZ basis sets (the aug-cc-pVTZ basis sets gave convergence problems). As shown in Fig. 4, one ZIF cluster and one CO₂ molecule were placed in different positions, then a set of configurations defined by moving the CO₂ molecule away from the ZIF cluster along a parallel direction. This time in the fitting we did not adjust the values of the ρ parameters but instead took these from the transferable potentials of Williams.²⁴ The final fitted parameters in Buckingham potential are given in Table 4.

Fig. 4 Two sets of CO₂ molecule and Zn(CN)₂ cluster configurations used in MP2 calculations.

Table 4 Fitted values of the parameters in Buckingham interaction (eqn (4)) between CO₂ molecules and the components of the ZIFs

Atom pair	B (kJ mol^−1)	ρ (Å)	C (kJ mol^−1 Å^−6)
CCO2–CZIF	106996.4	0.278	1487.38
CCO2–H	60 683	0.2793	713.33
CCO2–N	103 366	0.282	1473.04
CCO2–Zn	1 417 347	0.290	38015.38
O–CZIF	163 495	0.265	1394.10
O–H	57 298	0.266	613.98
O–N	157 948	0.269	1380.67
O–Zn	274 353	0.275	11435.99

---

3 Molecular dynamics simulations

Molecular dynamics (MD) simulations were performed using the DL_POLY code (version 4.06),³² running on the UK's Archer national high-performance computers. We used an ensemble with constant temperature and pressure controlled using the Nosé–Hoover thermostat and barostat, with a constraint to maintain 90° angles between the three axes of the atomic configuration. The simulations used time steps of 0.003 ps.

The simulation configurations consisted of slabs of the ZIF structures with three-dimensional periodic boundaries, with the space between the slabs filled with CO₂ molecules. The slabs were created always with the crystal [001] axis normal to slab surface, cutting the Zn–N bonds at a position designed to minimise the number of bonds to be cut. We always therefore align the normal to the slab with the sample Z axis, which will be denoted as such in the rest of this paper. The thickness of the slab was chosen to contain 10 cavities, and the lateral dimensions were chosen to be around twice the thickness. The thickness thus depended on the specific ZIF, and varied from 75 Å for ZIF-zni to 211 Å for ZIF-6. To make the system electrically neutral, we attached a hydrogen atom to each open nitrogen atom, as shown in Fig. 5.

Fig. 5 ZIF slabs used for CO₂ absorption simulation which have been run at temperature 200 °C, under pressure 25 bar. The Z axis, namely the axis normal to the plane of the slab, is vertical in each image. The imidizolate ligands are replaced by rods for clear illustration. The Zn atoms are denoted by tetrahedral shapes to show a better 3-D structure. The layer of pink atoms denotes the surface of the slab.

The simulations were performed under pressures and temperatures corresponding to the conditions that might be found in a typical industrial process. The plan with the simulations was that the pressure would force the CO₂ molecules to enter the ZIF slab. This would result in the volume occupied by the gas to be continuously decreasing, until almost all the molecules in that particular stage of the simulation were inside the slab. Typically the height of the gas column at the start of the simulation was around 800–1000 Å. After one simulation, which will have ended when most of the CO₂ molecules had entered the slab, we added a new column of CO₂ on top of the ZIF slabs once the CO₂ molecules get into the slab, as shown in Fig. 6. The left pane of the figure shows the start of the first stage simulation, and the right pane shows start of the second stage simulation with newly added CO₂ position above the slab and with CO₂ absorbed within the slab from the first stage simulation. We define one complete absorption process as the suite of simulations to the saturation state where no more molecules are easily absorbed. For one complete absorption process, several stage simulations are needed. The onset to saturation was most easily monitored through the rate of change of volume with time.

Fig. 6 One configuration of CO₂ on top of ZIF-3 slab at temperature 200 °C and pressure 25 bar. The imidazolate ligands are replaced by rods for clearer illustration. Zn atoms are denoted by tetrahedral shape to show a better 3D structure. A layer of pink atoms is the surface of the slab.

To perform the task of adding CO₂ molecules to the new DL_POLY configuration file and increase the dimension of the simulation box at each new stage during the absorption process we used a self-written Python script. During the simulation, the coordinates of the system shift along the Z axis when a large amount of CO₂ been pushed in, and this shift cannot be predicted, which means that the boundary of the simulation box might cut through the middle of the solid absorbent or the CO₂ gas. Due to this shift, any bonds that cross the periodic boundary along the normal axis would be elongated when the length of this axis be increased. Our python script was also used to divide the CO₂ gas by a given plane normal to the Z axis and reset this plane as the new periodic boundary.

For the data analysis, we wrote one Python script for extracting the output data from the file generated by DL_POLY that contains a set of data for each time step to be recorded (called STATIS), and another Python script to calculate the atom density profiles along the three axes. The visualisation software CrystalMaker³³ was used to produce snapshot images and movies to provide additional information.

In addition to the absorption simulation, reverse desorption simulations were also performed. These started with the configurations for which absorption had reached their equilibrium state, but the configuration was set to atmospheric pressure.

It is useful to note that the large volume changes that occur with out approach to study absorption and desorption are not really compatible with the way that the architecture of the MD software is optimised for parallel computing architectures. In particular it proved impossible to run desorption simulations with DL_POLY4, so for these simulations it was necessary to use the older DL_POLY Classic software, which uses a different parallelisation strategy more suitable for smaller samples. The computational overhead limited the extent to which it was possible to study the desorption process, so the results presented here are rather more limited than for the studies of absorption.

4 Simulations of the CO₂ absorption process

4.1 Absorption condition

Temperatures of 200 °C and 400 °C, and pressures of 25 bar and 30 bar, have been chosen as the absorption conditions use here following the information provided in two reports prepared for the International Energy Agency on oxy-combustion capture in power plants.^34,35 There are several temperatures and pressures used within the entire combustion process. Before combustion, high-pressure steam may be generated at pressures around 290 bar and temperatures around 600 °C. After combustion, the CO₂-rich gas may leave the heat recovery system at a temperature of 110 °C. The gas will then be compressed to 30 bar, half of which may be cooled to −55 °C for further compression for pipeline transmission, and the other half cooled to 35 °C to be recycled into the boiler. Considering the supercritical pressure of CO₂ is 75 bar rather than steam pressure, we chose to use pressure ranges from 25–30 bar for our simulations of the post-combustion absorption process. To represent the absorption under different temperatures, temperature ranges from 200 °C to 400 °C were chosen for our simulations. Based on the methods and results reported here, both the pressure and the temperature can be decreased or increased in future studies to represent the different operating conditions as desired.

4.2 Absorption at 25 bar and 200 °C

The absorption of CO₂ can be monitored through the system volume, since the simulation box will shrink along the Z axis as the gas molecule are pushed into the slab. We have performed a systematic study of CO₂ absorption in several ZIFs at a temperature of 200 °C and pressure of 25 bar. The variations with time of the volumes of several ZIFs during the absorption processes are shown in Fig. 7. In this figure, the Roman numerals against each curve denote the particular stage of the simulation. For example, in the case of ZIF-zni, the gas volume during stage I – the first stage starting with no CO₂ within the ZIF – decreased rapidly as CO₂ is pushed into the ZIF until almost all the CO₂ in the gas had been absorbed. Stage II then represents a re-start of the simulation from the end of stage I but with a new volume of gas added to the top. From the results we can see generally, and this is exemplified by the example of ZIF-zni, that the absorption process slows down as more CO₂ is pushed into the slab. ZIF-6 is an example that required more stages, with each stage taking longer than the previous stage; in this case stage I took around 500 ps, stage II took around 1000 ps, and later stages took even longer. When comparing results in Fig. 7 it should be noted that the same height and density of CO₂ molecules was added above every ZIF shown in Fig. 7. One exception to the general trends was ZIF-8, in that the volume change with time in ZIF-8 appears to have been almost linear, and within reasonable simulation time scales it was not possible to reach the equilibrium state. This appears to be due to the interaction between CO₂ molecules and imidazolate ligands; we will discuss this in more detail in a separate paper (manuscript in preparation).

Fig. 7 The change of volume during the CO₂ absorption process in several ZIFs at 200 °C under 25 bar from our molecular dynamic simulations. Each line labeled by the Roman numerals represents one simulation stage during which a volume of CO₂ has been pushed into the ZIF slab completely. Each stage run started from the ZIF configuration of the former stage run but with a new volume of CO₂ added above the ZIF slab.

The absorption results are summarised in Table 5, which characterises each ZIF in terms of the number of atoms in each simulation and the density of Zn atoms, and then comparing the number of CO₂ molecules absorbed in the simulation and the density of absorbed CO₂ molecules. For each case Table 5 also gives the time required to reach saturation in the simulation. The pore (or channel) density is inversely proportional to the density of Zn atoms, with a higher Zn density indicating a smaller pore structure and hence smaller pore surface areas. We might have expected to find a higher CO₂ gas capacity in the ZIFs with higher pore surface areas, but the simulation results in Table 5 show no simple link between the Zn density and CO₂ uptakes. The results implies that the uptake capacity is decided by more factors rather than just the Zn density.

Table 5 CO₂ absorption in ZIFs taken from the last outputs of the simulations. N_ZIF is the number of total atoms in ZIF slabs, N_CO2 is the total number of absorbed CO₂ molecules, D_Zn is the density of Zn atom in ZIF slabs, D_CO2 is the density of absorbed CO₂ molecules. The unit of the densities is mol ml⁻¹. The last column gives the time required for the CO₂ absorption process from the start to eventual saturation state

	N ZIF	N CO2	D Zn	D CO2	Time (ps)
ZIF-zni	28 008	1455	7.26	6.64	4859
ZIF-2	86 112	9781	5.35	10.64	6656
ZIF-3	80 896	11 464	4.30	10.69	1369
ZIF-4	80 896	9241	5.22	10.45	4029
ZIF-6	180 992	25 487	4.12	9.97	8205
ZIF-8	147 584	3553	4.06	2.30	21 085
ZIF-10	230 496	34 701	4.18	10.86	16 508

---

We remark that the density of the CO₂ gas outside the slab in our simulations after saturation with pressure of 25 bar and temperature of 200 °C is approximately 34 kg m⁻³. This is higher than both the experiment value of 29 kg m⁻³ at the same condition, which in turn is consistent with our own MD simulations of the gas state showing a density of 28.5 kg m⁻³ at the same conditions. We speculate that the increased density is associated with attractions due to the surfaces, although if this is the explanation it was not possible to test this because the whole sample size was not large enough to observe the expected gradient in the density in the direction between the surface of one slab and the closest surface of its periodic image.

4.3 CO₂ density profile

The CO₂ density profiles for the directions normal to the slab surfaces at different stage during the absorption process indicate that the gas molecules always accumulated in the surface layers at first and were then pushed further into the slabs by the accumulating gas molecules. We show this process for ZIF-2 in Fig. 8, In the first view, at 243 ps (first pane in the top row Fig. 8), the CO₂ molecules had accumulated in the layers near the surfaces (which are marked by the vertical dashed lines). With more gas molecules having been pushed into the surface, molecules already within the slap have been pushed into the deeper layers, as shown at 1524 ps (second pane in the top row and first in the bottom row in Fig. 8). This trend continues (third pane in the top row in Fig. 8) until the gas is saturated in the slab (the last two panes in Fig. 8). This trend is consistent with the rate of volume change of the gas over time (Fig. 7). From the start, the gas molecules only need to stay near the surface pores, and their capture takes place in short time. Subsequently they will be pushed further into the slab by the pressure from the gas stream, but this happens over a relatively longer period. Most of the other ZIFs we have studied, with the exception of ZIF-8, shown the same trends as ZIF-2, and their density profiles are given in the ESI.†

Fig. 8 Absorbed CO₂ density profile at 243 ps, 1524 ps, 3354 ps, 4605 ps, 6081 ps and 11 262 ps in ZIF-2, respectively from left to right and top to bottom. The Z axis is the axis normal to the slab. The two dashed lines denoted the edge of the slab and which define the average length between two layers.

The patterns of CO₂ absorbed in different ZIFs also reflect the structures of pores and channels. This is shown for example of ZIF-4 in Fig. 9. The observation of the pattern of pores reflects the high degree of filling that is being achieved in the simulation. The patterns for absorption in other ZIFs are given in the ESI.†

Fig. 9 CO₂ molecules in ZIF-4 slab after the simulations have reached the equilibrium state. The image shows a side view of CO₂ in the ZIF-4 slab without showing the atoms of ZIF-4. The connected blue tetrahedra are used to represent the pore structures of the ZIFs before absorption happens, and the large yellow spheres surrounded by the tetrahedra are used to represent the pore space.

4.4 Absorption at 30 bar

To further understand the absorption dynamics and to explore the CO₂ uptake capacity of ZIFs, the absorption process has been investigated by increasing the pressure from 25 bar to 30 bar. This pressure increase is a compromise between being large enough in order to see a difference, but not so large as to significantly increase the computational requirements. The starting configurations were taken from the last MD run at 25 bar in each case (saturation state under 25 bar).

The change in volume of the system (slab plus gas) in different ZIFs under increased pressure is illustrated in Fig. 10, where we show the effect of switching instantaneously from 25 bar to 30 bar. It can be seen that the volumes decreased rapidly under the higher pressures. This shrinking process will in part be associated with compression of the gas. However, there is also a higher degree of absorption, as shown by the results given in Table 6. They show that the absorption uptake capacity of most of the ZIFs (ZIF-zni, ZIF-2, ZIF-4, ZIF-6 and ZIF-10) under 30 bar is slightly higher than at 25 bar, by just under 5%. The increase in uptake is much larger in ZIF-3 and ZIF-8 slabs, of 8% for ZIF-3 and 16% for ZIF-8. This mechanism of this high increase can be revealed by the CO₂ density profiles within the ZIF. The density profile of absorbed CO₂ in the ZIF-3 slab under 30 bar is shown in Fig. 11. It shows that the density peaks for 30 bar are higher than the ones under 25 bar. This shows simply that more CO₂ molecules have been pushed in, and that they tend to cluster together in the same sites.

Fig. 10 Volume change of CO₂ absorption in ZIFs. The black lines represents the adsorption under 25 bar and the red lines represents the adsorption under 30 bar.

Table 6 CO₂ absorption in ZIFs under 25 bar and 30 bar. N_CO2 is the number of absorbed CO₂ molecules, D_CO2 is the density of the absorbed CO₂ molecules. It is given both in mol ml⁻¹ and kg m⁻³. P_increase is the increased percentage of the absorbed CO₂ molecules

	N CO2 (25 bar)	N CO2 (30 bar)	D CO2 (25 bar)	D CO2 (30 bar)	P increase
Units:			mol ml^−1/kg m^−3	mol ml^−1/kg m^−3
ZIF-zni	1448	1502	6.64/278.96	6.87/301.97	3.43%
ZIF-2	9739	10 013	10.64/468.12	10.90/479.55	2.42%
ZIF-3	11 458	12 403	10.69/470.21	11.52/506.82	7.84%
ZIF-4	9222	9462	10.45/459.75	10.68/469.87	2.22%
ZIF-6	25 371	26 561	9.97/438.63	10.41/547.99	4.47%
ZIF-8	3552	4066	2.30/101.19	2.63/115.7071	14.53%
ZIF-10	34 630	35 187	10.86/477.79	11.03/485.27	1.62%

---

Fig. 11 Absorbed CO₂ density profile in ZIF-3 slab. The black line represents the adsorption profile for 25 bar, the red lines represents the adsorption profile for 30 bar.

4.5 Interaction between imidazolate ligands and CO₂ molecules

We now look at the CO₂ absorption sites and consider the attraction between CO₂ and the imidizolate ligands. The N atoms in ZIFs are negatively charged,¹⁴ and the C atoms in CO₂ have positive charge. This leads to a Coulomb attraction between the centre of a CO₂ molecule and N_ZIF. The C_CO2–N_ZIF Coulomb attraction leads to the preferable adsorption sites that are illustrated in Fig. 12. It can be seen that most of the CO₂ molecules have become attached to the ligands, with the carbon atom in the CO₂ molecule being positioned close to the nitrogen atom in the imidazolate ligand. Inspection of Fig. 12 also shows that the CO₂ molecules tend to lie flat parallel to the plane of the imidazolite ligands, and evenly located on both sides.

Fig. 12 CO₂ adsorption sites in ZIF-zni, ZIF-2, ZIF-3, ZIF-4 an ZIF-6 pore and channels. H atoms are not shown for a better illustration. The grey atoms are Zn. Transparent light blue spheres are used to illustrate the pore structure. N and C atoms of the imidazolate ligands are represented by the wire bonds for a better illustration.

We have also made a animation of CO₂ in ZIF-zni at the saturation state, using CrystalMaker, which is provided within the ESI.† The animation shows that the CO₂ molecules in saturated ZIF-zni move back and forth with the translational and rotational vibrations of the ligands.

5 CO₂ desorption from ZIFs

5.1 Introduction

Having shown that most of the ZIFs (ZIF-8 being the key exception) will absorb CO₂ from a gas with a modest pressure until some saturation state, we have next to ask the question as to whether the process can be reversed. For applications of ZIFs for carbon capture, being able to recover the CO₂ from the ZIF to enable subsequent storage is essential. Traditionally there are two ways to encourage desorption, namely increasing the temperature or decreasing the pressure. We decrease the pressure from 25 bar to 1 bar whilst maintaining the same temperature of 200 °C, starting with the configuration corresponding to the equilibrium state of saturation absorption. In addition to testing decreasing pressure at the same temperature, we also have tried to simulate the desorption in ZIF-10 by increasing the temperature to 400 °C.

The way to handle the desorption simulation is basically a reverse process of absorption simulation. Generally, the simulations have been performed in separate stages similar to absorption simulations. When the gas was released to a certain amount, typically an expansion in the direction of the normal to the slab by a factor of about 4 or 5, the gas molecules were removed from the simulation manually and the volume adjusted accordingly. This new configuration was then the starting point for the subsequent simulation.

The desorption process was initially simulated using the same MD code DL_POLY used for the absorption simulations, but we found that this code gave some problems associated with its use of high-performance computing architectures. We found that the same problems did not arise using a previous version, DL_POLY classic (version Classic 1.9), but this has the disadvantage of not being optimised for large systems in the same way that Version 4 has been. Thus we used version 4.06.1.1 for the initial desorption simulations, and then we switched to DL_POLY classic for further desorption runs. This strategy, however, only allowed us to perform short period simulations of the desorption process.

5.2 Desorption dynamics and mechanisms

As with the absorption studies, in our simulations on the desorption process we monitored the volume as representing the extent of gas release, as shown in Fig. 13. Again, different stages are indicated by Roman numerals. The first stage, denoted by I, corresponds to the start of the desorption process, in which we started with the sample of saturated systems obtained at pressure of 25 bar and temperature of 200 °C and then reduced the pressure to 1 bar. The volume increased as CO₂ was released from the slab. Subsequent stages (II, III etc.) follow after the manual removal of molecules from the gas as described above. It can be seen that the volume of the systems expand rapidly at I stage. The desorption process then slows down as the volume expands, which is exactly the opposite to the absorption process. Take the desorption in ZIF-3 for instance; starting from the same volume, the first run (I) took 40 ps, but the second run (II) took approximately 160 ps to expand to the same volume. The desorption processes in ZIF-2, ZIF4 and ZIF-10 have shown the similar trend as ZIF-3. Because of this we have only been able to study the desorption processes at their early stage.

Fig. 13 CO₂ desorption dynamics in ZIFs at 200 °C under 1 bar. The configuration of the first run (I run) is taken from the last output files of absorption under 25 bar.

The decreased percentage of absorbed CO₂ molecules are shown in Table 7. The Table shows that big portions of the gas have been released in ZIF-2, ZIF-3 and ZIF-4 already. It takes longer time to release the gas in ZIF-2. However, their trends in Fig. 13 indicates that the systems are still far away from the equilibrium state which means a larger amount of CO₂ can still be released from these ZIFs, including ZIF-10.

Table 7 CO₂ desorption in ZIFs under 1 bar. N_CO2 is the number of absorbed CO₂ molecules, time (ps) is the time elapsed since the pressure changed into 1 bar from 25 bar. Decrease percentage is the (released CO₂ amount from the slab at 1 bar)/(absorbed CO₂ amount at 25 bar)

	N CO2 (25 bar)	N CO2 (1 bar)	Time (ps)	Decreased percentage
ZIF-2	9781	8710	402	16.47%
ZIF-3	11 464	8547	555	25.44%
ZIF-4	9241	6588	903	28.71%
ZIF-10	34 701	32 978	462	5%

---

We have also performed a simulation of the desorption process at the higher temperature of 400 °C, with results shown in the last image in Fig. 13. It took approximately 66 ps at 200 °C to expand to the same volume, but only 42 ps at 400 °C. The number of CO₂ molecules that came out of the slab further showed that the gas molecules are released quicker at higher temperature. While only 684 CO₂ molecules have been released at 200 °C after 66 ps, in the 42 ps to give the same volume the slab released three times the number of molecules, namely 2294 molecules. This means that the absorbed CO₂ molecules can be released significantly faster by heating up the system.

6 Conclusion

In this paper we have studied the absorption of CO₂ in zeolitic imidazolate framework (ZIF) structures using the molecular dynamics simulation method. In order to do this work, we have built on our prior simulation study of the ZIF structures for which we had developed a significant part of the force field,¹⁴ and for the present paper we augmented this with a new force field for CO₂ designed to work in both condensed and fluid phases together with new models for the interactions between CO₂ and ZIF. The new model for CO₂ has been validated with calculations of crystal structure and phonon dispersion curves, and of the liquid–gas coexistence behaviour.

Based on these model force fields, a set of systematic MD simulations on CO₂ absorption process in ZIF slabs have been presented. The simulations on CO₂ absorption in ZIFs have been performed with environments close to those found in industrial power plants. The simulations have produced a large quantity of data, some of which have been reported in this paper but with many results presented in the ESI.† The simulations have given insights into the process of absorption, showing both similarities and differences between the behaviour of different ZIF systems. It is found that it is usually possible to achieve a high filling of the porous structures, as monitored both by calculations of density profiles and visualisations. It has been found that the preferable absorption sites for CO₂ molecules are close to the imidazolate ligands due to the C_CO2–N_ZIF Coulomb attraction.

In addition to the absorption simulation, the desorption processes have also been simulated. The desorption simulation results have shown that the CO₂ molecules can be taken out by decreasing the pressure or increasing the temperature.

All these results point to the real possibility of using porous structures like ZIFs for capture of CO₂ in an industrial environment as part of a strategy to reduce the release of CO₂ into the Earth's atmosphere. In a following paper we will discuss some of the details of the processes at the atomic scale, but here our headline results concern the potential of such materials for this application, but with the caveat that different members of the general family of materials may show very different absorption capacity.

In future work it would be useful to work with gas streams containing other molecules, such as N₂ and H₂O. This will require an extension of our force field model, ideally using the approaches discussed in this paper.

Acknowledgements

This research utilised two high-performance computing facilities. Development of the force field was carried out using Queen Mary's MidPlus computational facilities, supported by QMUL Research-IT and funded by EPSRC grant EP/K000128/1. The molecular dynamics simulations were carried out using the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk), with access made available through our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202). MG and CY were supported by both the China Scholarship Council and Queen Mary University of London. AM was supported by a European Union Marie Sklodowska-Curie fellowship.

References

1. D. Venkataraman, S. Lee, J. Zhang and J. S. Moore, Nature, 1994, 371, 591–593.
2. O. M. Yaghi and H. Li, J. Am. Chem. Soc., 1995, 117, 10401–10402.
3. J. L. Rowsell and O. M. Yaghi, Chem. Eng. J., 2004, 73, 3–14.
4. K. Sumida, D. L. Rogow, J. A. Mason, T. M. McDonald, E. D. Bloch, Z. R. Herm, T.-H. Bae and J. R. Long, Chem. Rev., 2012, 112, 724–781.
5. K. S. Park, Z. Ni, A. P. Cote, J. Y. Choi, R. Huang, F. J. Uribe-Romo, H. K. Chae, M. O'Keeffe and O. M. Yaghi, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 10186–10191.
6. R. Banerjee, A. Phan, B. Wang, C. Knobler, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Science, 2008, 319, 939–943.
7. A. Phan, C. J. Doonan, F. J. Uribe-Romo, C. B. Knobler, M. O'Keeffe and O. M. Yaghi, Acc. Chem. Res., 2010, 43, 58–67.
8. Y. Li, F. Liang, H. Bux, W. Yang and J. rgen Caro, J. Membr. Sci., 2010, 354, 48–54.
9. Q. Song, S. K. Nataraj, M. V. Roussenova, J. C. Tan, D. J. Hughes, W. Li, P. Bourgoin, M. A. Alam, A. K. Cheetham, S. A. Al-Muhtaseb and E. Sivaniah, Energy Environ. Sci., 2012, 5, 8359–8369.
10. S.-S. Han, S.-H. Choi and W. A. Goddard III, J. Phys. Chem. C, 2011, 115, 3507–3512.
11. U. P. N. Tran, K. K. A. Le and N. T. S. Phan, ACS Catal., 2011, 1, 120–127.
12. A. Phan, C. J. Doonan, F. J. Uribe-Romo, C. B. Knobler, M. O'Keeffe and O. M. Yaghi, Acc. Chem. Res., 2010, 43, 58–67.
13. L. T. L. NGuyen, K. K. A. Le and N. T. S. Phan, Chin. J. Catal., 2012, 33, 688–696.
14. M. Gao, A. J. Misquitta, L. H. N. Rimmer and M. T. Dove, Dalton Trans., 2016, 45, 4289–4302.
15. S. S.-Y. Chui, S. M.-F. Lo, J. P. H. Charmant, A. G. Orpen and I. D. Williams, Science, 1999, 283, 1148–1150.
16. E. Aprà, E. J. Bylaska, D. J. Dean, A. Fortunelli, F. Gao, P. S. Krstić, J. C. Wells and T. L. Windus, Comput. Mater. Sci., 2003, 28, 209–221.
17. M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. J. Van Dam, D. Wang, J. Nieplocha, E. Apra, T. L. Windus and W. de Jong, Comput. Phys. Commun., 2010, 181, 1477–1489.
18. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
19. A. J. Stone, The Theory of Intermolecular Forces, Oxford University Press, 2nd edn, 2012, pp. 122–140.
20. A. J. Misquitta and A. J. Stone, CamCASP 5.9, 2016, http://www-stone.ch.cam.ac.uk/documentation/camcasp/users_guide.pdf.
21. P. J. Winn, G. G. Ferenczy and C. A. Reynolds, J. Phys. Chem. A, 1997, 101, 5437–5445.
22. G. G. Ferenczy, P. J. Winn and C. A. Reynolds, J. Phys. Chem. A, 1997, 101, 5446–5455.
23. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
24. D. E. Williams, J. Comput. Chem., 2001, 22, 1154–1156.
25. B. L. Gastaldi and L. Scaramuzza, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1980, 36, 2751–2753.
26. D. Bolmatov, D. Zav'yalov, M. Gao and M. Zhernenkov, J. Phys. Chem. Lett., 2014, 5, 2785–2790.
27. J. J. Potoff and J. I. Siepmann, AIChE J., 2001, 47, 1676–1682.
28. A. Simon and K. Peters, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1980, 36, 2750–2751.
29. G. Dolling, P. Martel, L. Piseri, P. Martel and B. Powell, Bull. Am. Phys. Soc., 1972, 17, 291.
30. E. W. Lemmon, M. O. McLinden and D. G. Friend, Thermophysical Properties of Fluid Systems in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, ed. P. J. Linstrom and W. G. Mallard, National Institute of Standards and Technology, Gaithersburg MD, 20899,  DOI:10.18434/T4D303, (retrieved July 21, 2017).
31. M. G. Martin, Mol. Simul., 2013, 39, 1212–1222.
32. I. T. Todorov, W. Smith, K. Trachenko and M. T. Dove, J. Mater. Chem., 2006, 16, 1911–1918.
33. D. C. Palmer, Z. Kristallogr. - Cryst. Mater., 2015, 230, 559–572.
34. IEA Greenhouse Gas R&D Programme (IEA GHG), Oxy-combustion processes for CO₂ capture from power plant, Report number 2005/9, July 2005. http://www.ieaghg.org/docs/General_Docs/Reports/Report%202005-9%20oxycombustion.pdf, retrieved July 2017.
35. IEA Greenhouse Gas R&D Programme (IEA GHG), CO₂ capture ready plants, Report number 2007/4, May 2007, https://www.iea.org/publications/freepublications/publication/CO2_Capture_Ready_Plants.pdf, retrieved July 2017.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7me00034k

‡ Present address: China Spallation Neutron Source, Institute of High Energy Physics, Chinese Academy of Sciences, 1 Zhongziyuan Road, Dongguan 523803, People's Republic of China, and, Dongguan Institute of Neutron Science, 1 Zhongziyuan Road, Dongguan, 523808, People's Republic of China; E-mail: E-mail: gaomin@ihep.ac.cn

§ Present address: Institute of Natural Sciences, and, Department of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, People's Republic of China.

¶ Present address: Fachbereich Physik, Erwin-Schrödinger-Strasse, 67663 Kaiserslautern, Germany.

---