Phonons and Adsorption-Induced Deformations in ZIFs: Is It Really a Gate Opening?

: We report a microscopic model of the phonon-adsorption correlations in ﬂ exible metal − organic framework materials. We analyze the mechanism of the gate opening deformation using the notion of coupled phonon- and adsorption-induced structural transformation. Using the ZIF-8 structure as an example, we perform an analysis of transformation-related, low-frequency phonon modes of the framework. On the basis of structure-related quantities such as pore limiting diameter, void fraction, and adsorption uptake, we determine the conditions which lead to the gate opening transformation in ZIF-8. Energetic landscape of the deformation process is analyzed using grand thermodynamic potential of adsorption. We generalize our conclusions to other ﬂ exible ZIF structures with the same topology.


■ INTRODUCTION
Metal−organic frameworks (MOFs) 1 are crystalline, hybrid (organic−inorganic) nanoporous structures with tunable architecture. They are potentially applicable in many fields of chemistry and chemical engineering, such as gas capture, separation, catalysis, sensing, and many others. 2 A subgroup of MOFs, flexible MOFs, are frameworks which undergo a variety of structural deformations induced by different external stimuli such as pressure, temperature, electric field, and adsorption. 3 To explain the mechanism of adsorption-induced transformations, the phonon-driven deformation approaches have been explored; 4−6 nevertheless, despite many numerical and experimental works devoted to flexible MOFs, a comprehensive, generalized model of the microscopic mechanism of MOFs' deformations is still missing, mainly due to a variety of origins of the observed phenomena (rotation or buckling of the linkers, 7−10 sublattice displacement, 11 change of the angle between linkers, 12 and so forth).
ZIF-8, a sodalite (SOD) topology MOF built of Zn cations interconnected by imidazole linkers (Figure 1, center), undergoes a transformation usually referred to as "gate opening", involving a rotation of linkers which enlarges (or diminishes) the pore aperture. 10 This naturally impacts the adsorption-related properties, as the pore aperture (the gate) defines the ability of guest molecules to diffuse through the pores of the framework. In ZIF-8, the gate opening was previously reported to occur under high hydrostatic pressure 13−15 and was related to the adsorption of N 2 . 10 However, the discrimination as to whether the step-on isotherm of nitrogen adsorption in ZIF-8 is caused by gate opening (linker reorientation) or by reorganization of the Article pubs.acs.org/JPCC adsorbate within the pore into a more densely packed configuration without adsorbent transformation was never achieved. 9,10 Most of the current theoretical models assume that the step present on the adsorption isotherm is associated with deformation of the host framework, defined by the frozen, low energy (frequency) phonons corresponding to linkerrotation-like vibrations. [4][5][6]16 Here, we analyze these phonons and discuss their role in the transformation occurring in ZIFtype MOFs.
In this work, we applied phonon formalism, density functional theory (DFT), and Grand Canonical Monte Carlo (GCMC) adsorption calculations to revisit the mechanism of the structural transformations in ZIF-type MOFs with SOD topology. We performed DFT calculations of low-frequency phonons and analyzed their potential contribution to the gate opening mechanism. For the first time, analysis of phonondriven deformations of ZIFs was combined with quantities describing adsorption process such as pore limiting diameter (PLD) and void fraction (VF), explicit adsorption simulations, and grand thermodynamic potential of adsorption. This accounts for both phenomena in one model.

■ COMPUTATIONAL METHODS AND MODELS
Preparation of Phonon-Deformed Structures. In this work, we analyze structures which are derived from a DFToptimized unit cell, where atoms are displaced according to deformation vectors u(i) related to particular phonon modes. These vectors are derived from polarization vectors e α (α is Cartesian coordinate) being eigenvectors of a 3D eigenvalue problem 17 where ω(q) is the frequency of a mode at given q-point and D(q) is the dynamical matrix. The deformation vector represents the instantaneous position of atom corresponding to a particular phonon where m i is the mass of ith atom, Q is the amplitude of mode (scale of the eigenvector used to create deformed structure; we refer to it as deformation parameter,) and e iqR is the plane wave which introduces the spatial dependence. Since we consider only phonons from the Γ-point of the Brillouin zone (for conventional, cubic cell q = 0), then the last term in eq 2 is equal to unity, and the deformation vector can be rewritten as In this paper, the eigenvectors e α are extracted from VASP output file (OUTCAR) and deformed structures are prepared with our in-house python script. 18 The deformation parameter Q was varied from −1.5 to 1.5 [in units of Å·amu 1/2 ] with step of 0.1. The value Q = 0 corresponds to an undistorted (optimized) ZIF-8 structure, whereas for the limiting value Q = 1.5 we observe full opening of the gates in a previously reported gate opening mode 4−6 (with linkers in a four-member ring arranged in parallel to each other).
DFT Calculations. All phonon calculations have been performed on optimized structures for which the internal energy was found to be minimal. The structure optimization was performed using the VASP package and PAW method. 19−21 As previously reported, 22 long-range interactions were taken into account using Grimme D3(BJ) correction 23 to PBE density functional. 24 The convergence criteria were set to 10 −6 eV for electronic degrees of freedom and to 10 −4 eV/Å for ionic. Energy cutoff was set to 900 eV. Only one (Γ) kpoint was considered due to large conventional unit cell of ZIF structures. Optimization was performed without symmetry constraint; the final symmetries of the optimized frameworks are listed in the Table S1. Frequencies are calculated using the finite displacement method with six displacements for each atom (±Δx, ±Δy, ±Δz, where each displacement was set to 0.01 Å), also considering only Γ-point (sufficiently large unit cell).
GCMC Calculations. GCMC simulations of the adsorption were performed using RASPA software. 25 Structures of the frameworks were taken from DFT calculations. To describe interactions in our system we used Lennard-Jones potential (6-12) parametrized according to UFF+ force field 26 for framework atoms and TraPPE 27 for N 2 molecules. UFF+ is based on UFF force field where ε parameters of LJ potential are divided by the factor of 1.69. Cross parameters were obtained with Lorentz−Berthelot mixing rules. Charges of framework atoms are calculated with eQeq 28 (method implemented in RASPA), as suggested in the benchmarking article by Sladekova et al. 29 Charges for N 2 molecule were taken from TraPPE force field. Spherical cutoff for LJ potential was set to 14 Å at which specific interactions were truncated (not shifted), and all interactions beyond this distance were calculated within mean-field approach (i.e., tail-correction was applied).
PLD and VF were calculated using algorithms available in the ZEO++ software. 30 All atomic structures and accessible surfaces were visualized in VESTA 31 and iRASPA 32 software, respectively. Absolute change of PLD and VF is defined as where X is considered parameter. Accessible Surface. Accessible surfaces were calculated using iRASPA software 32 with nitrogen as a probe molecule. This surface is an isosurface where energy of interactions between probe and framework is equal to 0 (r = σ LJ ). Practically, it means that the volume inside the pore corresponds only to the positions of the center of mass of adsorbed molecules. Hence, even a very thin channel which connects the neighboring pores indicates the possibility of gas diffusion in the framework.
Potential of Adsorption. Grand thermodynamic potential 33 for adsorbed phases is given as ads 0 (5) where N(p,Q,T) is the adsorption isotherm calculated for a unit cell deformed with deformation parameter Q and temperature T.

■ RESULTS AND DISCUSSION
Analysis of the DFT-derived phonon modes and their polarization vectors leads us to the conclusion that the lowfrequency range of the vibrational spectra in ZIF-8 is totally dominated by modes related to linker rotations. The modes are divided into two groups, symmetric and asymmetric (for the We focus our analysis on two important modes occurring at 64.6 and 111.7 cm −1 . The mode at frequency 64.6 cm −1 (we called it ν go1 mode, Figure 1, top) was previously considered to be the most important in the gate-opening transformation 4−6 and for Q = 1.5 this mode changes the angle between linkers in the four-member ring by around 24°, leading to an essentially parallel orientation of the linkers. However, to the best of our knowledge no explicit relation between the deformation caused by the mode and its impact on the experimental adsorption isotherm was ever established. To get better insight into this aspect, we analyzed two structural parameters sensitive to the phonon deformations and, at the same time, affecting adsorption properties of a given structure: PLD, the diameter of the smallest pore aperture, and VF, the fraction of an empty space in volume of porous material. Both parameters were calculated for stable configuration of ZIF-8 (with geometry optimized at T = 0 K), and for hypothetical structures deformed along the deformation vectors u α of few modes from low-frequency regime. The effect of deformation for selected modes and Q = ±1.5 is shown on Figure 1. We follow the convention adopted in the literature, 4,6 and call the deformations leading to the increase of the access to the pore a positive deformation. Figure 2 shows the dependence of PLD and VF on the deformation parameter Q. Out of all modes listed in Table 1, the four modes impacting PLD and VF in the most significant way are shown here. Surprisingly, deformation related to the ν go1 mode causes a decrease of both parameters by 30% and 24% respectively. In consequence, the pore access is blocked, and the pore volume is decreased; consequently, the deformation rather should be called a "gate closing" deformation. The mode at frequency 103.8 cm −1 is doubly degenerated, however, as the deformations related to this mode exceeds by far harmonic region and nonzero deformation removes degeneracy, the corresponding VF and PLD are slightly different.
To visualize this transformation better, Figure 3 shows the surface accessible for adsorption in undistorted and distorted ZIF-8 structures. In undistorted structure, neighboring pores are connected by narrow apertures (formed by six Zn atoms connected by imidazole linkers, six-member ring), and therefore the diffusion of the gas through the system is possible. In the structure distorted according to ν go1 mode, these apertures disappear. It is a consequence of the blocking of six-member ring gates by hydrogen atoms from methyl groups attached to rotated linkers. As deformation progresses, for Q ≥ 1 the PLD (Figure 2) decreases indicating that the gate is closing. A similar process is observed for the mode with frequency 76.95 cm −1 at even smaller deformation parameter Q ≥ 0.5.
The mode with the frequency of 111.7 cm −1 (we call it ν go2 mode, see Figure 1) also involves linker rotation but has totally the opposite effect on PLD and VF parameters as compared to ν go1 mode. It is symmetric with respect to the negative and positive deformation and causes a constant, monotonic increase of PLD and VF as the deformation progresses (by 16% and 7%, respectively, for Q = 1.5, Figure 2). Consequently, both negative and positive deformations are in fact the same (graphical explanation provided in Figure S2). The increase of the gate diameter and subsequent variation of geometry of surface accessible for adsorption in the vicinity of the pore apertures ( Figure 3, right panel) increases the probability of the guest molecules to diffuse between the neighboring pores. Therefore, ν go2 mode appears as a precursor of adsorption-induced gate opening deformation.
To analyze the explicit correlation between the deformation and adsorption, we performed simulations of nitrogen adsorption (at T = 77 K in undistorted ZIF-8 structure and in structures distorted according to ν go1 and ν go2 modes. Several simulation unit cells were prepared for deformation   The Journal of Physical Chemistry C pubs.acs.org/JPCC Article parameter varying from 0 to 1.5. Figure 4 shows isotherms of adsorption for the outermost cases (Q = 1.5) compared to that calculated for undistorted structure (the isotherms for all intermediate structures (0 < Q < 1.5) for both modes are presented in the Figure S1) and an experimental isotherm. 9 For structures deformed along ν go1 mode, the adsorbed amount is smaller than in undistorted structure in the entire range of studied pressures. This observation suggests that ν go1type of deformation (i) decreases the strength of adsorption sites (at low pressures) and (ii) decreases the pore volume (at high pressure), as compared to the undistorted structure. In the structure deformed by ν go2 mode, the adsorbed amount is larger than in the undistorted structure (in the whole pressure range) although smaller than the one measured in the experiment. 9,10 Experimental adsorption isotherm (Figure 4, black) shows two steps: at the pressures of 0.005 and 0.009 bar. Experimentally, the first step was attributed to the gateopening mechanism and the second to an expansion of the ZIF-8 framework (after structural deformation) upon adsorption. 9 The discrepancy between experiment and approach presented in this work can be related to several factors. First, in our calculations we use a ZIF-8 unit cell which, prior to phonon and adsorption calculations, was optimized at T = 0 K. Since ZIF-8 is a flexible material with positive volumetric thermal expansion coefficient, 35 we expect that at finite temperatures its unit cell and hence VF enlarges, causing an increase of the maximal gas uptake. Furthermore, in finite temperatures there are more types of vibrations occurring simultaneously, which can lead to a modification of the total capacity of the framework. Figure 5 shows the change of grand thermodynamic potential of adsorption ΔΩ ads as a function of structure deformation (0 < Q < 1.5) according to ν go1 and ν go2 modes. ΔΩ ads provides insight into the system energy variation upon adsorption-induced deformation and was calculated from nitrogen adsorption isotherms ( Figure S1; for theoretical details, see SI). For ν go1 -like deformation, ΔΩ ads rises monotonically as the deformation increases, reaching the value of almost 70 kJ/mol (of ZIF-8 conventional unit cell) for Q = 1.5 (fully open gate structure), hence this process requires the energy to be supplied to the system. On the contrary, in the case of ν go2 -like deformation ΔΩ ads progressively decreases, reaching the value of −35 kJ/mol for fully open gate conformation. As the energy required to deform ZIF-8 to a gate-opened phase, evaluated using ab initio molecular dynamics, is equal to approximately 15 kJ/mol, 36 we conclude that ν go2 -deformed structure is energetically stable in the presence of guest molecules adsorbed in the pore of ZIF-8, unlike the ν go1 -deformed structure. In other words, the ν go1like deformation may not be induced by adsorption.
To check whether chemical modification of ZIF-8 structure affects the mechanism of gate-opening transformation and related adsorption properties, 4,14,15,37 we extended our analysis to other ZIF materials which share the same sodalite topology and differ either in (i) functional group attached to imidazole linker (ZIF-90, 38 SALEM-2, 39 ZIF-8(Amino), 40 and ZIF-65 41 ) or (ii) metal ion (BIF-3(Li), 42 BIF-3(Cu), 42 ZIF-67, 43 and CdIF-1 44 ). The adopted selection protocol and a detailed description of frameworks is given in SI. Figure 6 (top) gives the values of frequencies of ν go1 and ν go2 modes in modified frameworks. As both modes exist in all studied structures, we conclude that they are characteristic for sodalite topology of the framework. The changes of PLD and VF in all distorted structures with respect to undistorted ones are given in Figure  6 (middle and bottom, respectively). For all frameworks with modified linkers, ν go1 and ν go2 related deformations result in a   These results suggest that ZIF-type structures with SOD-type topology are more prone to undergo gate opening transformation initiated by ν go1 and or ν go2 lattice vibrations when the imidazole linker of the initial ZIF-8 structure is chemically modified. Two cation-modified frameworks (ZIF-67 and CdIF-1) show similar response to the ν go1 and ν go2 deformations as ZIF-8 (PLD is smaller or larger, respectively). For all systems deformed along ν go1 mode, the VF decreases, whereas for those deformed along ν go2 increases by several percent. An exception from this rule is observed for BIF-type materials of smaller (by 17%) unit cell volume than ZIF-8, resulting from shorter imidazole-B bonds (1.6 Å, to be compared to imidazole-Zn bond of 2.0 Å in ZIF-8).

■ CONCLUSIONS
In conclusion, our work shed new light on the question of the origin of the structural transformations in ZIF-type MOFs with sodalite topology. For the first time, analysis of phonon-driven deformations of ZIFs was combined with simulations of adsorption, to account for both phenomena at the same time.
We have shown that low-frequency vibration mode at ν go1 = 64.6 cm −1 , which was previously associated with deformation causing an opening of the gates between pores actually causes gates' closing and disruption of the connections between neighboring pores. In consequence, the pore volume accessible for adsorption decreases in disaccord with the experimental observation of a rapid increase of the adsorbed amount at structure transformation. Moreover, in the low-frequency regime of ZIF-8 phonon spectrum we have identified another vibration mode at ν go2 = 111.7 cm −1 , consisting of a rotation of the imidazole linkers, which enlarges a six-member ring forming the pore aperture. The deformation of the structure according to this mode progressively increases both the diameter of the gate between the pores (gate opening deformation), and the fraction of structure volume accessible for adsorption. In consequence, the adsorption amount in the deforming structure is larger than in undeformed one which stays in agreement with the experimental results. On the basis of analysis of grand thermodynamic potential of adsorption, we prove that this type of the deformation is energetically favored. We also showed that the linker-rotation-type deformation of the framework is typical for imidazolate-based structures that show sodalite-topology but differ in the functional group attached to linker and metal cation. Since in solids there are many types of motions occurring simultaneously, in further studies the combination of two or more vibrational modes should also be considered as a source of static deformations of flexible structures. Overall, the presented study confirms that the phonon analysis coupled with simulations of adsorption constitutes a useful tool to interpret adsorption-induced deformations and should be widely considered as a source of valuable microscopic information about the systems, complementary to experimental explorations.
Details on simulation methodology used in this work and description of the linker-and cation-modifications to ZIF-8 structure; animations of linker rotation modes; animation of the accessible surface construction (PDF)