Aurélie U Ortiz

Anne Boutin

Alain H Fuchs

François-Xavier Coudert

Investigating the Pressure-Induced Amorphization of Zeolitic Imidazolate Framework ZIF-8: Mechanical Instability due to Shear Mode Softening

Keywords: metal-organic frameworks, shear instability, mechanical properties, molecular simulation, elastic constants, structural transition, nanoporous materials = shear instability ZIF-8 + pressure C44

We provide the first molecular dynamics study of the mechanical instability that is the cause of pressure-induced amorphization of zeolitic imidazolate framework ZIF-8. By measuring the elastic constants of ZIF-8 up to the amorphization pressure, we show that the crystal-to-amorphous transition is triggered by the mechanical instability of ZIF-8 under compression, due to shear mode softening of the material. No similar softening was observed under temperature increase, explaining the absence of temperature-induced amorphization in ZIF-8. We also demonstrate the large impact of the presence of adsorbate in the pores on the mechanical stability and compressibility of the framework, increasing its shear stability. This first molecular dynamics study of ZIF mechanical properties under variations of pressure, temperature and pore filling opens the way to a more comprehensive understanding of their mechanical stability, structural transitions and amorphization.

Zeolitic imidazolate frameworks (ZIF) are a subclass of metal-organic frameworks (MOF) [START_REF] Horike | Soft Porous Crystals[END_REF][START_REF] Farha | Rational Design, Synthesis, Purification and Activation of Metal-Organic Framework Materials[END_REF][START_REF] Férey | Hybrid Porous Solids: Past, Present, Future[END_REF][START_REF] Tranchemontagne | Secondary Building Units, Nets and Bonding in the Chemistry of Metal-Organic Frameworks[END_REF] that feature imidazolate linkers bridging metal centers to form three-dimensional porous crystalline solids isomorphous to zeolitic frameworks. [START_REF] Park | Exceptional Chemical And Thermal Stability Of Zeolitic Imidazolate Frameworks[END_REF][START_REF] Zhang | Pore Surface Tailored SOD-Type Metal-Organic Zeolites[END_REF][START_REF] Phan | Structure, And Carbon Dioxide Capture Properties Of Zeolitic Imidazolate Frameworks[END_REF][START_REF] Zhang | Metal Azolate Frameworks: From Crystal Engineering To Functional Materials[END_REF] ZIFs have recently gained considerable attention for their potential applications because they inherit desirable qualities from both worlds: the tunable porosity, structural flexibility and the functionalization of the internal surface of the MOFs, as well as the thermal, mechanical and chemical stability of the zeolites. Moreover, the similarities between the metalimidazolate four-fold coordination chemistry and the corner-sharing SiO4 tetrahedra from which zeolites are built mean that many ZIF topologies can potentially be synthesized. Indeed, over 100 different ZIF structures have been reported so far. Theoretical [START_REF] Lewis | Zeolitic Imidazole Frameworks: Structural And Energetics Trends Compared With Their Zeolite Analogues[END_REF][START_REF] Baburin | Enumeration Of Not-Yet-Synthesized Zeolitic Zinc Imidazolate Mof Networks: A Topological And Dft Approach[END_REF] as well as experimental [START_REF] Hughes | Thermochemistry of Zeolitic Imidazolate Frameworks of Varying Porosity[END_REF] investigations of the relative stabilities of these structures indicate that they have a relatively low enthalpy of formation, and that many more polymorphs should be accessible under mild synthesis conditions. The ZIF structures reported in the literature so far demonstrate large gamut of attractive structural and physicochemical properties, with great potential for applications in fluid separation, [START_REF] Li | Zeolitic Imidazolate Frameworks For Kinetic Separation Of Propane And Propene[END_REF][START_REF] Gücüyener | Ethane/Ethene Separation Turned on Its Head: Selective Ethane Adsorption on the Metal-Organic Framework ZIF-7 through a Gate-Opening Mechanism[END_REF][START_REF] Bux | Oriented Zeolitic Imidazolate Framework-8 Membrane with Sharp H2/C3H8 Molecular Sieve Separation[END_REF][START_REF] Bae | A High-Performance Gas-Separation Membrane Containing Submicrometer-Sized Metal-Organic Framework Crystals[END_REF] CO2 capture, [START_REF] Banerjee | High-Throughput Synthesis of Zeolitic Imidazolate Frameworks and Application to CO2 Capture[END_REF] sensing, [START_REF] Lu | Metal-Organic Frameworks As Sensors: A ZIF-8 Based Fabry-Pérot Device As A Selective Sensor For Chemical Vapors And Gases[END_REF] encapsulation and controlled delivery. [START_REF] Liedana | CAF@ZIF-8: One-Step Encapsulation Of Caffeine In MOF[END_REF] Given their larger number of potential applications, issues of thermal and mechanical stability of flexible MOFs in general, [START_REF] Ortiz | Anisotropic Elastic Properties of Flexible Metal-Organic Frameworks: How Soft are Soft Porous Crystals?[END_REF][START_REF] Ortiz | Metal-Organic Frameworks With Wine-Rack Motif: What Determines Their Flexibility And Elastic Properties?[END_REF][START_REF] De Toni | How Can a Hydrophobic MOF be Water-Unstable? Insight into the Hydration Mechanism of IRMOFs[END_REF] and ZIFs in particular, [START_REF] Tan | Exceptionally Low Shear Modulus in a Prototypical Imidazole-Based Metal-Organic Framework[END_REF] are of high importance. A number of recent studies have addressed the question of ZIF stability at high temperature, showing temperature-induced amorphization of some ZIFs (ZIF-1, ZIF-3, ZIF-4), [START_REF] Bennett | Structure and Properties of an Amorphous Metal-Organic Framework[END_REF] while others do not amorphize at high temperature (including the prototypical ZIF-8). [START_REF] Bennett | Thermal Amorphization of Zeolitic Imidazolate Frameworks[END_REF] On the other hand, all ZIF frameworks studied so far have shown pressure-induced amorphization at modest, industrially accessible pressure. ZIF-8 [START_REF] Park | Exceptional Chemical And Thermal Stability Of Zeolitic Imidazolate Frameworks[END_REF] (Fig. 1) undergoes irreversible amorphization under nonhydrostatic mechanical compression as low as 0.3 GPa, [START_REF] Chapman | Pressure-Induced Amorphization and Porosity Modification in a Metal-Organic Framework[END_REF] as well as under ball-milling. [START_REF] Cao | Amorphization Of The Prototypical Zeolitic Imidazolate Framework ZIF-8 By Ball-Milling[END_REF] The pressure-induced amorphization of ZIF-4 [START_REF] Park | Exceptional Chemical And Thermal Stability Of Zeolitic Imidazolate Frameworks[END_REF] (Fig. 1) was also studied: Bennett et al. [START_REF] Bennett | Reversible Pressure-Induced Amorphization Of A Zeolitic Imidazolate Framework (ZIF-4)[END_REF] showed that amorphization occurs at various values of hydrostatic pressure depending on the presence or absence of solvent molecules in the pores of the material, as well as the nature of the hydrostatic medium. The evacuated ZIF-4 demonstrated amorphization at very low pressure (0.35 GPa), while the presence of DMF molecules in its pores could shift amorphization to higher pressure and lead to the existence of an intermediate monoclinic crystalline phase (ZIF-4-I). While all the experimental studies of pressure-induced amorphization cited above have taken great care in identifying the structure of the various amorphous ZIF phases, there is a lack of experimental information on the nature of the crystal-to-amorphous transition, the reasons behind its occurrence and its apparent universality under mechanical pressure. This is in good part because the amorphization happens at such low pressure that little data is available on the stressed ZIFs. Furthermore, a recent study by Tan et al. [START_REF] Tan | Exceptionally Low Shear Modulus in a Prototypical Imidazole-Based Metal-Organic Framework[END_REF] measured the shear modulus of ZIF-8, which they find to be quite low, at around 1 GPa. Tan et al. also suggested that the pressure-induced amorphization of ZIF-8 may be linked to this low shear modulus. However, the pressure-induced amorphization is observed at pressure lower than 1 GPa, and in addition a direct link between these two properties has yet to be demonstrated. In this work, we report the evolution of the mechanical properties of ZIF-8, as a function of both external stress and adsorbate loading, in order to elucidate the nature of the pressure-induced crystal-to-amorphous transition.

ZIF-8 (SOD)

ZIF-4 (cag)

N N N N Zn + N N N N Zn +
Figure 2 presents the evolution of the unit cell parameter and elastic constants of ZIF-8 as a function of mechanical pressure, at 300 K, for an isotropic compression corresponding to hydrostatic experimental conditions with a non-penetrating hydrostatic fluid. Both the average lattice parameters and the elastic constants were obtained from constant-stress (N, σ, T) molecular dynamics (MD) simulations by analysis of the fluctuations of unit cell vectors [START_REF] Ray | Fluctuations And Thermodynamic Properties Of Anisotropic Solids[END_REF][START_REF] Parrinello | Strain Fluctuations And Elastic Constants[END_REF] (the so-called strain-fluctuation method), where the elastic stiffness tensor Cijkl is obtained by:

! " B # $ % 𝐶 '(") *+ = 〈𝜖 '( 𝜖 ") 〉 -〈𝜖 '( 〉〈𝜖 ") 〉 (1)
with ε the unit cell strain and V the unit cell volume . We described the fully flexibly ZIF-8 framework using the force field of Zhang et al, which was calibrated to reproduce the ZIF-8 structure, bulk modulus and sorption-induced structural transformation [START_REF] Zhang | Sorption-Induced Structural Transition of Zeolitic Imidazolate Framework-8: A Hybrid Molecular Simulation Study[END_REF] (see Supporting Information for simulation details). First, direct observation of the MD trajectories shows that the ZIF-8 framework is stable at pressures lower than 0.4 GPa, but unstable at higher pressure. In the range of pressures from -0.5 GPa to 0.35 GPa (where negative pressures indicate tension and positive pressures indicate compression), ZIF-8

shows linear elasticity with a bulk modulus of E = 8 GPa, in good agreement with the values determined experimentally (from Brillouin scattering: 7.7 GPa; [START_REF] Tan | Exceptionally Low Shear Modulus in a Prototypical Imidazole-Based Metal-Organic Framework[END_REF] from pressure-induced lattice variations:

6.5 GPa 25 ). At P ≥ 0.4 GPa, however, we observe a rapid collapse of the framework, demonstrating its [START_REF] Chapman | Pressure-Induced Amorphization and Porosity Modification in a Metal-Organic Framework[END_REF] To find the fundamental cause of this instability we look at the evolution of the elastic constants of ZIF-8 as a function of applied pressure (Fig. 2, bottom panel). Because the ZIF-8 framework is cubic, its mechanical properties are fully described by only 3 elastic constants: C11, C12 and C44, which are respectively the modulus for axial compression, dilation on compression, and shear. First, we see that the values obtained for the unstressed crystal (P = 0) using this methodology are in good agreement with the experimental values obtained by Brillouin scattering: the agreement is very good for C11 and C12 (this work: 11.3 and 7.6 GPa; ref. 22: 9.5 and 6.9 GPa), and reasonable for C44 (2.7 GPa vs. 0.97 GPa) considering how challenging is the accurate prediction of absolute values of elastic constants. More important is the pressure-dependence of the elastic constants, where we can see that C44 diminishes drastically under compression of the crystal, a phenomenon known as shear-mode softening. This decrease leads to elastic instability of the ZIF-8 framework, since the Born stability conditions for a cubic crystal under hydrostatic pressure P can be written as: [START_REF] Morris | The Internal Stability Of An Elastic Solid[END_REF] C11 + 2C12 + P ≥ 0 ; C11 -C12 ≥ 2P ; C44 ≥ P

(2)

At P = 0.4 GPa the third condition, which indicates stability under shear deformation, breaks down and the ZIF-8 framework becomes mechanical unstable. We thus conclude that the amorphization of ZIF-8 under pressure is triggered by elastic instability, itself due to a pressure-induced softening of the (already relatively soft) shear deformation mode of the material. While it is the first time such a mechanism is demonstrated for amorphization of a MOF, it is known to be the underlying cause of pressureinduced phase transitions in some dense inorganic solids (e.g. the halite-CsCl-type transition of MgO, or the stishovite-CaCl2-type transition in SiO2 32 ) as well as the pressure-induced amorphization of coesite. [START_REF] Dean | Pressure-Induced Amorphization In Crystalline Silica: Soft Phonon Modes And Shear Instabilities In Coesite[END_REF] In their recent study of ZIF-4 amorphization, Bennett et al. [START_REF] Bennett | Reversible Pressure-Induced Amorphization Of A Zeolitic Imidazolate Framework (ZIF-4)[END_REF] have demonstrated that both the onset of amorphization and the bulk modulus of ZIF-4 depend highly on the presence or absence of solvent molecules in its pores. In particular, they showed that the evacuated framework was much more flexible.

Moreover, it was shown in some soft porous crystals, such as the MIL-53 family, that adsorptioninduced stress can lead to structural transitions of large amplitude. [START_REF] Triguero | Understanding Adsorption-Induced Structural Transitions In Metal-Organic Frameworks: From The Unit Cell To The Crystal[END_REF] These considerations are of critical importance for practical applications in real conditions. We have thus investigated the elastic constants and mechanical stability of ZIF-8 with methane adsorbed inside its pores. Methane was chosen as a generic guest molecule for this study because of its apolar nature and the absence of specific guest-ZIF interactions, in order to better understand the generic effect of pore volume filling on mechanical properties. Table I reports the evolution of ZIF-8 elastic constants for CH4 loadings between 0 and 18 molecules per unit cell (which corresponds to methane uptake at 100 bar and 300 K; see isotherm in Fig. S1).

While the effect of pore filling is modest on the C11 and C12 elastic constants, it drastically increases the shear modulus of ZIF-8, with a variation of C44 from 2.7 GPa in the evacuated material to 4.0 GPa with 18 methane molecules per unit cell. This increased resistance to shear is directly linked to the higher density of the guest-filled material, for which a shear motion implies a larger number of repulsive guest-guest and guest-host close contacts. red: empty ZIF-8 at 300 K; blue: empty ZIF-8 at 500 K; orange: ZIF-8 with 18 CH4 at 300 K.

Moreover, MD simulations of ZIF-8 at this loading and under hydrostatic compression have revealed its increased mechanical stability compared to the empty ZIF-8: the material is stable under compression up to 0.65 GPa, while it collapses at higher pressure. The pressure dependence of the elastic constants is shown in Fig. 3 in the form of the stability criteria for compression (λ = C11 -C12 -2P) and shear (μ = C44 -P), both of which must be positive for the system to be mechanically stable. Fig. 3 demonstrates that pressure-induced amorphization in methane-loaded ZIF-8 is, like the evacuated material, triggered by a mechanical instability due to shear softening, but that the onset of amorphization is delayed due to the larger shear modulus of the unstressed material. We also plot on Fig. 3 the evolution of elastic constants at 77 K and 500 K. No temperature-induced softening is observed: C44 is almost unchanged, while C11 actually increases. The absence of elastic softening at higher temperature is in good agreement with the fact that no amorphization of ZIF-8 upon heating has been reported in the literature.

Preliminary results on the elastic constants and mechanical stability of the ZIF-4 framework under hydrostatic pressure (Fig. S2) point to the same mechanism for its pressure-induced amorphization. A softening of the shearing modes (three of them because the ZIF-4 framework is orthorhombic) is also found to be associated with mechanical instability under compression, though the value of transition pressure obtained by this study does not fall in line with the experimental results and further work needs to be done to refine the force field used for the simulation to account for the differences in coordination bonds and angles between ZIF-8 (for which the force field was optimized) and ZIF-4. Nonetheless, the fact that the very different ZIF-8 and ZIF-4 frameworks both undergo pressure-induced amorphization due to shear instability hints that the phenomenon may be generic. We thus conclude that the porous ZIF phases, though they retain zeolite-like framework topologies, do not possess the mechanical stability of zeolites because of the weaker Zn-imidazolate coordinative bonds, compared to the strong Si-O bond.

Moreover, their very high porosity gives rise to low resistance to shear and pressure-induced shear softening, which in turn is the cause for their mechanical instability at low pressure.

Finally it should be noted that, while the molecular dynamics simulations reported in this work are helpful to understand how the mechanical instability of ZIFs develop depending on thermodynamic conditions, they unfortunately cannot yield direct information on the nature of the amorphous phase obtained and do not explain why the amorphous phase is preferred to a hypothetical denser crystalline phase. All simulations of ZIF-8 at high pressure, in which the framework collapses, lead to the occurrence of a denser porous crystalline phase of Zn(mim)2, with reduced pore volume compared to ZIF-8.

This can be attributed to the limitations of the force field used: by its bonded nature, it does not allow reconstructive transitions to take place in the system. DFT-based molecular dynamics simulations would allow a direct simulation of the amorphization process, but the computational cost of such simulations on a large-scale supercell of ZIF-8 is not within our reach.

In conclusion, we have determined the mechanism of the pressure-induced amorphization of ZIF-8 as a pressure-induced mechanical instability due to shear mode softening. The shear modulus (C44) of the ZIF-8 cubic framework, which was shown by earlier work to be quite low, actually decreases under compression of the framework. We find that temperature changes do not induce instability of the framework, consistently with the absence of temperature-induced amorphization of ZIF-8. We show that the presence of guest molecules in the highly porous ZIF-8 greatly increases its resistance to shear, requiring a higher pressure to amorphize. In essence, we highlight how constant-stress molecular dynamics can be used to investigate pressure-and temperature-induced amorphization and structural transitions.

ASSOCIATED CONTENT

Supporting Information. Molecular simulation details, values of elastic constants, methane adsorption isotherms. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 1 .

 1 Figure 1. Imidazolate linker, structure and topology of ZIF-4 and ZIF-8. Yellow spheres represent the

Figure 2 .

 2 Figure 2. Unit cell parameter a and elastic constants C11, C12 and C44 of ZIF-8 as a function of hydro-

  under pressure. The value of this critical pressure is, again, in good agreement with the experimental data of Chapman et al. who observed amorphization beyond 0.34 GPa under hydrostatic conditions.

Figure 3 .

 3 Figure 3. Evolution of the stability criteria for compression (upper panel, λ = C11 -C12 -2P) and shear

Table I . Elastic constants of the ZIF-8 framework at 300 K for various values CH4 loading, at 300 K and 0 GPa.

 I

	CH4 / unit	C11	C12	C44
	cell	(GPa)	(GPa)	(GPa)
	0	11.21	7.49	2.75
	3	11.26	7.52	2.89
	6	11.29	7.55	2.96
	11	11.33	7.58	3.18
	18	11.38	7.66	4.00

AUTHOR INFORMATION

Corresponding Author * fx.coudert@chimie-paristech.fr

ACKNOWLEDGMENT

Funding from the Agence Nationale de la Recherche under project "SOFT-CRYSTAB" (ANR-2010-BLAN-0822) is acknowledged. This work was performed using HPC resources from GENCI-IDRIS (projects 100254 and 100167), the NAMD software, [START_REF] Phillips | Scalable Molecular Dynamics with NAMD[END_REF] which was developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana-Champaign.


Volker Haigis
email: volker.haigis@ens.fr

Xavier Coudert

Rodolphe Vuilleumier

Anne Boutin

Investigation of structure and dynamics of the hydrated metal-organic framework MIL-53(Cr) using first-principles molecular dynamics

The hydration behavior of metal-organic frameworks (MOFs) is of interest both from a practical and from a fundamental point of view: it is linked, on the one hand, to the hydrothermal stability (or instability) of the nanoporous material, which might limit its use in technological applications. On the other hand, it sheds light on the behavior of water in a strongly confined environment.

Here, we use first-principles molecular dynamics (MD) to investigate two hydrated phases of the flexible MOF MIL-53(Cr), which adopts a narrow-or a large-pore form, depending on the water loading. Structure and dynamics of the two phases are thoroughly analyzed and compared, with a focus on the hydroxyl group of MIL-53(Cr) and the water molecules in the nanopores. Furthermore, the behavior of the confined water is compared to that of bulk water. Whereas in the narrow-pore form, water is adsorbed at specific crystalline sites, it shows a more disordered, bulk-like structure in the large-pore form. However, reorientation dynamics of water molecules in the latter is considerably slowed down with respect to bulk water, which highlights the confinement effect of the nanoporous framework.

Introduction

Metal-organic frameworks (MOFs) are a topical class of nanoporous materials displaying a large range of crystal structures and host-guest properties, due to a combination of tunable porosity, by choice of metal centres and linker length, and functionalisation of the internal surface of the material. They have gained a lot of attention in the last decade and have been proposed for use in applications such as adsorptive storage, gas separation, catalysis and sensing. However, even though a rapidly increasing number of structures have been synthesized, there are far less detailed characterizations of water adsorption in MOFs in the literature [1][2][3][4][5][6][7][8][9] than adsorption studies on other gases of strategic interest, such as carbon dioxide, methane and hydrogen. Nonetheless, the behavior of MOFs in presence of water vapor, their stability under humidity and the impact of water adsorption on their other physical and chemical properties (adsorptive separation performance, catalysis, etc.) are of primary importance if one envisions applications at the industrial level.

Our interest in the present work is to characterize the structure, vibrational and orientational dynamics of water confined in the archetypical material MIL-53(Cr). This material is formed of unidimensional chains of corner-sharing CrO 4 (µ 2 -OH) 2 octahedra, linked by 1,4-benzenedicarboxylate (BDC) ligands to form linear diamond-shaped channels large enough to accommodate small guest molecules 10 (Fig. 1). This structure may oscillate between two phases: a large-pore (lp) structure and a narrow-pore (np) one. This phenomenon, which has been termed "breathing", can be induced by adsorption or desorption of guest molecules including water 11,12 , variations in temperature 13 or mechanical stress 14 . Full studies of adsorption and desorption of water, among other polar vapors, in MIL-53(Cr) were published by Bourrelly et al. 15 and Devautour-Vinot et al. 16 including experimental sorption characterization, calorimetric measurements and quantum chemistry calculations of the preferential arrangements of the molecules within the pores. These quantum chemistry calculations, performed at the density functional theory 17,18 (DFT) level, were geometry optimizations leading to water structures of lowest energy for each phase of the host (narrow and large pore). A later experimental study of MIL-53(Cr) immersed in liquid water used high resolution powder X-ray diffraction to solve the structure of the pores of both the large-pore and narrow-pore phases, fully filled with water 19 .

In addition, recent molecular simulation studies have used force field-based molecular dynamics and grand canonical Monte-Carlo simulations to study the structure, diffusion and adsorption thermodynamics of water in the lp and np phases of MIL-53(Cr) 20,21 . The two works differ in their specifics: Salles et al. 20 focused on classical (Newtonian) dynamics of the water and framework, while Paesani 21 used normal-mode path integral molecular dynamics in order to include quan-tum dynamical effects. However, both studies described the interactions in the system at the classical level, with the use of a parameterized nonpolarizable force field. Furthermore, Salles et al. 20 describe the water molecule as a rigid body. Both approximation are quite drastic, and may limit the precision with which water-MOF and water-water interactions are described. Indeed, a recent report by Cirera et al. 22 demonstrated that the effects of electronic polarization significantly impact the structure of water in the pores of MIL-53(Cr), especially at large loading of the unit cell (fully hydrated pores). Very recently, Chen et al. 23 published an ab-initio study on the structural response of the Sc-bearing MIL-53(Sc) to temperature changes and CO 2 adsorption.

Here, we present a study of the structure and dynamics of the fully hydrated MIL-53(Cr) np and lp phases, using firstprinciples molecular dynamics. We used Car-Parrinello molecular dynamics, meaning that the interactions in the system are treated at the quantum mechanical level, in the DFT approach, while the molecular dynamics samples the phase space of the system and allows the calculation of statistical averages and dynamical properties (and not only minimal energy configurations, as in quantum chemistry energy minimizations). We have shown earlier that this approach allows a good description of the structure, dynamics, electronic properties and reactivity of molecular fluids confined in porous spaces. In particular, we showed how it could describe the behavior of water confined in a family of zeolites 24,25 and explain the limited hydrothermal stability of MOFs from the IR-MOF family 9 . Compared to simulations based on empirical force fields, the first-principles molecular dynamics approach has the advantage that it does not depend on empirical parameters and hence is transferable to very different chemical environments, e.g. situations where chemical reactions occur. Moreover, in the present case of nano-confined water, the validity of classical force fields has to be checked carefully since they are parameterized for bulk water in general and might not perform well in surface-dominated systems. With the firstprinciples approach, we aim at providing a reliable reference and at validating previous results obtained from force-field simulations. Also, it paves the way for future studies of chemical reactivity in these systems. However, the high computational costs of this approach impose limits on the system size and simulation times, and phenomena such as gas diffusion in nanoporous materials may still lie beyond the capabilities of first-principles simulations.

Computational methods

The hydrated phases of MIL-53(Cr) were studied by means of DFT-based Car-Parrinello molecular dynamics simulations as implemented in the CPMD package 26,[START_REF] Marx | Modern Methods and Algorithms of Quantum Chemistry[END_REF] . The simulation cell represents one unit cell of the empty (non-hydrated) MIL-53(Cr), containing 4 [Cr(OH)-(O 2 C-C 6 H 4 -CO 2 )], and the hydrated narrow-and large-pore phases were modeled by adding 4 and 24 water molecules, respectively, to the empty framework. This corresponds to the amount of water found experimentally in the two phases 19 . The simulations were performed in the canonical (NV T ) ensemble, with experimental unit cell parameters taken from Guillou et al. 19 . The temperature of the ions was maintained at 350 K by a Berendsen thermostat with a coupling time constant of 0.5 ps. It was shown that time constants larger than 0.1 ps yield reliable structural and dynamical properties of water [START_REF] Berendsen | [END_REF] . Another Berendsen thermostat was used to constrain the kinetic energy of the fictitious degrees of freedom with a coupling time constant of 0.25 ps. The mass associated with the latter was chosen to be 200 a.u., and a time step of 0.05 fs was used for the integration of the equations of motion. The interactions between electrons and ion cores were described by Troullier-Martins pseudopoten-tials 29 , and the PBE exchange-correlation functional 30 was used. Dispersion interactions, which are known to be poorly represented by the semilocal PBE functional, were included via the scheme proposed by Grimme 31 . A plane-wave cut-off of 120 Ry was found to be sufficient to obtain converged ionic forces. All simulations were started from the structures of hydrated MIL-53(Cr) as determined by x-ray diffraction (narrow-pore and large-pore form, Guillou et al. 19 ). Since the hydrogen atoms are not visible by x-ray diffraction, we added them "by hand" in reasonable positions. For compu-tational convenience, the heavier deuterium isotope was used throughout the study. Its larger mass slows down the hydrogen dynamics and allows using a larger time step, without affect-ing the chemical properties of the system. The limited size of the simulation cell is due to the high computational demands of firstprinciples simulations and might affect properties such as the kinetics of phase transitions or diffusion coefficients, and these phenomena are therefore not discussed in this study. However, the finite-size effect on structural, vibrational, and rotational properties is expected to be small.

In addition to the simulations of MIL-53(Cr), we generated a trajectory of 13 ps for bulk water, represented by 64 H 2 O molecules in a periodically repeated simulation cell and with a density of 1 g/cm 3 . The settings for this simulation was the same as for MIL-53(Cr), except for a larger mass of 400 a.u. of the fictitious degrees of freedom, which were not thermostat-ted, and a larger time step of 0.073 fs.

Results and discussion

Structure of the narrow-pore and large-pore phases

The hydrated narrow-pore phase with four water molecules in the simulation box, named MIL-53(Cr)(1H 2 O per Cr), was equilibrated for 5 ps and subsequently simulated for another 30 ps. During the whole simulation, the oxygen atoms of the water molecules formed hydrogen bonds with the hy-droxyl groups of the MIL-53 framework. No hopping of water molecules to neighboring binding sites was observed within the time span of the MD run. This suggests that the wa-ter inside the narrow pores may be considered as part of the crystal structure rather than a fluid phase. To corroborate this finding, we averaged, over the whole trajectory, the positions of all atoms, except the hydrogen atoms belonging to water molecules, and checked the resulting structure for symmetry. It was found that the average structure has indeed space group C2/c, the same as the non-hydrated MIL-53(Cr). The wa-ter hydrogen atoms were excluded from this analysis because they rotate freely about the dipole axis of the water molecules, and hence their average positions are not meaningful. The average structure is shown in Fig. 1.

The simulation of the hydrated large-pore phase with 24 wa-ter molecules in the simulation box, named MIL-53(Cr) (6H 2 O per Cr), was started from the structure determined by x-ray diffraction 19 . The system was equilibrated for 5 ps and sub-sequently simulated for another 25 ps. For the Rietveld refinement of the X-ray data, the authors imposed symmetry constraints on the positions of the water molecules. How-ever, in our simulations, no symmetry was assumed, and all atoms were completely free to move within the orthorhombic simulation box. This enabled us to check whether the water molecules within the pores really occupy well-defined crystalline sites or if they form a more disordered structure, possibly resembling bulk water. Indeed, a disordered hydration structure was found in this case, which will be disussed in section 3.4.1. Regarding the average structure of the framework itself, a tilting of the organic ligands with respect to the plane defined by neighboring inorganic chains is observed, which was not seen in the x-ray diffraction. Because this tilting is one of two symmetrically equivalent tilts, it is thus absent from the experimental crystallographic structure, which is an aver-age (non-tilted) structure of the two conformations.

Behavior of hydroxyl groups

We then turned our attention to the dynamics of the hydroxyl group of the Cr(µ 2 -OH) inorganic chain of the MIL-53(Cr) framework. While its oxygen atom is rather constrained and has a low amplitude of movement, with an average mean square displacement of about 0.2 ˚A with respect to its average position (see Table 1), the position of the hydrogen atom (and thus the orientation of the µ 2 -OH group) vary more widely during the course of our simulations. In both the lp and np phases, the average position of the H atom lies in the CrOCr plane, but exhibits large wagging motions around this symmetric position. We quantified these by plotting the distributions of the angles between the µ 2 -OH group and the CrOCr plane during the molecular dynamics (Fig. 2). It can be seen that in both structures, the wagging motion of the µ 2 -OH group can reach angles of up to 40°, with a slightly larger amplitude for the narrow-pore phase than in the large-pore phase (respective full half width of the distribution: 33°and 27°; obtained from a Gaussian fit). We attribute this difference to larger steric hindrance in the lp phase, due to the high density of water molecules inside the porous channels. A similar conclusion can be reached on the basis of the mean square displacement of the H atom from its mean position, with values of 0.314 Å and 0.305 Å in the np and lp phases. This large amplitude of the wagging motion is consistent with the fact that the proton position could not be solved in the experimental crystallographic structures. Remarkably, the Gaussian-like nature of the µ 2 -OH wagging angle reveals that this vibration mode is harmonic in nature, even with such a large amplitude. species np phase lp phase

O µ 2 -OH 0.200 Å 0.207 Å H µ 2 -OH 0.314 Å 0.305 Å O w 0.811 Å 0.484 Å Table 1
Mean square displacement for selected species with respect to their average position, averaged over equivalent atoms. For the lp phase, only the hydroxyl groups and the water molecules which remain hydrogen-bonded to each other during the whole simulation are taken into account.

Water molecules bound to µ 2 -OH groups

After characterizing the µ 2 -OH group of the MIL-53(Cr) material, we now focus on the water molecules hydrogen-bonded to these groups. As discussed in section 3.1, the situation here is very different between the narrow-pore and large-pore phases. In the np phase, there is a 1:1 stoichiometry between water and µ 2 -OH, with each water molecule bound to an hy- droxyl group, with O-O distances and O-H-O angles typical of hydrogen bonds. No hopping of the water was observed during the 30 ps of dynamics. However, while the water molecules do not diffuse, they show large orientational freedom: the molecules rotate rather freely in a half-space, while keeping intact the hydroxyl-water hydrogen bond. As a part of this rotational diffusion, we observe during the dynamics many instances where a water's H atom comes in close vicinity of a neighboring carboxylate group. In the closest case, we observed O water -O carboxylate distances down to 2.5 Å, and H water -O carboxylate distances down to 1.6 Å. However, such events cannot be construed as proper hydrogen bonds on account of their short lifetime (typically less than 200 fs). A similar phenomenon is known to take place for water confined in siliceous zeolites, and is not due to specific water-material interactions but merely to the narrow confinement of the water molecule 25 .

In the large-pore phase, the situation is quite different. Out of the four µ 2 -OH groups in the unit cell, three remain hydrogen-bonded to the same three water molecules throughout the molecular simulation. For the fourth hydroxyl, however, we do observe an exchange of the hydrogen-bonded water molecule twice during the dynamics. This can be seen in Fig. 3 seen at t = 4 ps and 7 ps. Even after 7 ps, the plot still reveals some large fluctuations of the O-O distance, with peaks up to 5 Å. This weaker hydrogen bonding is also reflected in the histogram of O µ 2 -OH -O water distances (Fig. 3, lower panel), which shows a long tail at large distances for the lp phase, in contrast with the narrow-pore case. This distribution of O µ 2 -OH -O water distances can be related to the vibrational spectrum of the µ 2 -OH group, which was obtained from the MD trajectories by Fourier-transforming the velocity autocorrelation function of µ 2 -OH bond vibrations (Fig. 4). The spectra of the narrow-and the large-pore form both represent an average over all four hydroxyl groups in the simulation cell. Since in the simulations, hy-drogen was represented by deuterium isotopes, the wavenum-bers of µ 2 -OH vibrations lie below the usual values around 3400 cm -1 . While the narrow-pore phase exhibits a well-defined peak around 2400 cm -1 , the spectrum of the large-pore phase is considerably broadened. This finding is in line with the more disordered hydration structure found in the lat-ter (see section 3.4.1) and is a direct consequence of the wider distribution of µ 2 -OH -O water distances in the large-pore phase (Fig. 3). The weaker hydrogen bonding resulting from the pro-nounced tail of the O µ 2 -OH -O water distances leads to a large-wavenumber shoulder in the µ 2 -OH stretching vibrations be-tween 2500 cm -1 and 2700 cm -1 , which is not observed in the narrow-pore phase.

At this point, we would like to comment on a point discussed by Cirera et al. 22 who used classical force fields to investigate the hydrogen bonding between water molecules and the µ 2 -OH groups in hydrated MIL-53(Cr). They found that a polarizable force field leads to a substantially different structure than a non-polarizable model. In particular, with a polarizable force filed, Cirera et al. 22 obtained a pronounced peak in the radial distribution function for water oxygen atoms and the H atoms of µ 2 -OH groups just below 2 ˚A (Fig. 4 in their study). This peak was weakened -and even absent at high water loadings -when the non-polarizable force field was used. The authors concluded that polarizability is essential in order to capture correctly the water-framework interaction. Our DFT simulations provide a suitable reference to check this statement, since the electronic degress of free-dom are treated explicitly and hence polarizability is automat-ically taken into account. Fig. 5 shows the radial distribu-tion function for O µ 2 -OH -O water in the narrow-and large-pore phase. Indeed, we obatained a strong peak around 1.8 ˚A in both cases, confirming the importance of polarizability to describe water-framework interactions accurately. However, we also found significant differences between the first-principles and the force-field simulations: at high water loading of MIL-53(Cr), the radial ditribution function goes to zero at distances around 3 ˚A (Fig. 5) in our simulations, whereas it is almost flat above 2.5 ˚A with the polarizable force field of Cirera et al. 22 .

Finally, we look at the environment of the µ 2 -OH-bound water molecules in the large-pore phase, and their hydrogen bonding pattern to neighboring molecules. Like in the case of the np phase, they can get close to the carboxylate groups of the framework, without forming H bonds. However, because of the presence of more water molecules in the larger pores, they do form H bonds with other water molecules. These hydrogen bond networks are naturally dynamic and live for few picoseconds, but two specific arrangements are repeatedly observed: a 2A1D ("A": accepting, "D": donating) configura- tion, in which the µ 2 -OH-bound water molecule is H-bond donor to one other molecule, and H-bond acceptor for the µ 2 - OH group and one H 2 O; a 1A2D, where the µ 2 -OH-bound water molecule is H-bond donor for two neighboring water molecules. Snapshots of these typical arrangements are depicted in Fig. 6. This is to be contrasted with the typical "2A2D" configuration in bulk water, where a water molecule accepts two and donates two H bonds.

3.4 "Bulk-like" water in the large-pore phase 3.4.1 Structure Fig. 7 shows the radial distribution functions for oxygen atoms which are part of water molecules (O water , upper panel) and for O -H pairs (lower panel). For comparison, bulk water is also shown. The broad peaks indicate a liquid-like structure. The first and second peak in the O water -O water distribution are located at 2.8 Å and 4.5 Å, respectively, in fair agreement with experimental data on bulk water at ambient conditions 32 as well as with the bulk water simulation results. Regarding the radial distribution of O water -H pairs, three peaks were found in the simulation: the first one at 1 Å corresponds to intramolecular bonding, the second one around 1.8 Å falls into the range of typical hydrogen bond lengths, and the third peak at 3.3 Å comprises non-bonded O water -H pairs. Note that for this analysis, all hydrogen atoms were taken into account, including those which are part of the MIL-53 framework. While the overall shape of the distribution functions resembles that of bulk water, the water inside the framework exhibits higher peaks and therefore appears to be more strongly structured.

While radial distribution functions reveal important aspects of liquid structure, they represent an isotropic average and are not unique, in the sense that very different atomic structures can lead to the same radial distribution functions. Additional information is contained in the distribution of angles formed by triplets of atoms. Fig. 8 shows the distribution of angles defined by three neighboring O water , and for reference, the bulk water distribution is plotted. The distribution forms a broad peak around 105 • , close to the tetrahedral angle of 109.5 • , with a minor feature at 50 • . The distance cutoff for neighboring O water was set to 3.5 Å, which corresponds to the first minimum of the respective radial distribution function, see Fig. 7. The angle distribution indicates that the O water form a network of more or less distorted tetrahedra, as expected for bulk water and also seen in the bulk water simulations. However, water in the MIL-53(Cr) framework shows a narrower distribution than bulk water, which corroborates the slightly overstructured nature of the confined liquid. The peak at small angles (in both the and the confined water phase) results from a combined excluded-volume and entropy effect: O 1 w -O 2 w -O 3 w angles smaller than 45°generate too short O 1 w -O 3 w distances. On the other hand, small angles are entropically favorable because they increase the available configuration space for other neighboring O water . In a compromise, water configurations accumulate around the smallest possible angle, which leads to the peak at 50°. In Fig. 8, we also plot the distribution of the angles between the intramolecular O water -H bond and the hydrogen bond formed by the same H. A cutoff of 1.25 Å and 2.5 Å was used for intramolecular and hydrogen bonds, respectively, as determined from the minima of the radial distribution functions, see Fig. 7. A broad peak is found which exhibits a maximum at 165 • . The distributions of confined and bulk water are of very similar shape.

Water reorientation dynamics

Given the structural similarities of water confined in the large-pore phase and bulk water, one might prematurely conclude that the behavior of the confined water is already close to the bulk limit. However, this limit is characterized not only by structure but also by dynamical properties. We therefore investigated the reorientation dynamics of water inside the MIL-53(Cr) framework by means of the correlation function

C(t) = u(0) • u(t) (1)
where u(t) is the normalized bisector of the H-O-H angle of a water molecule at time t, and the average . . . is taken over all water molecules and different time origins. In order to obtain information on possible anisotropic dynamics, we also calculated the component-wise correlation functions

C α (t) = u α (0)u α (t) (2)
with α ∈ {x, y, z}. They sum up to the total C(t) at every instant of time. The results for the large-and narrow-pore phases are shown in Fig. 9. In the large-pore phase, the total correlation functions decays much more slowly than in bulk water (dashed line). For the time interval 1 ps < t < 10 ps, it can be fitted by an exponential C(t) = A exp(-t/τ) with a time constant τ = 32 ps. For bulk water, a value of τ = 13 ps is obtained, i.e. we find a ratio τ lp /τ bulk = 2.5. The component-wise correlation functions reveal anisotropic behavior: the values of C α (t = 0) indicate a preferential orientation of the H-O-H bisector along z and y, i.e. along the diamond-shaped channels of the framework and in the direction joining opposite µ 2 -OH groups (see Fig. 1 for the labelling of the axes). Furthermore, the z and y components decay more slowly than the x component and are thus responsable for the retarded decay of the total correlation function with respect to the one of bulk water. We note that at the time scale of several nanoseconds, which is not accessible by DFT-based MD, Salles et al. 20 found relaxation times up to 2 ns, using an empirical force field and water different from ours.

In the narrow-pore phase, on the other hand, C(t) decays very rapidly for short times t < 1 ps. At larger times, the highly anisotropic nature of the reorientation dynamics prevents us from reducing it to a single time constant by an exponential fit. The component-wise correlation functions reveal that again, the water orientation along y and z is strongly preferred. While C x and C z decrease rapidly, C y remains constant after an initial decay. This reflects the stability of the water ori-entation, pointing away from the µ 2 -OH groups of the frame-work, and we anticipate that the constant C y will eventually dominate the total C(t) at large t. Hence a significantly slower relaxation is expected for large time scales.

In our approach, the atomic nuclei follow a classical trajectory on a potential energy surface determined by quantum interactions. In general, the quantum nature of the nuclei, in particular of the protons in water molecules, can strongly affect the dynamical properties of a system. However, in the present case, the effect of quantum dynamics seems to be largely reduced by the confinement of water in nanopores: Paesani 21 found that relaxation times for the reorientation of water molecules in MIL-53(Cr), derived from classical and quantum dynamics, agree to within a factor 1.3 for most water loadings, and always to within a factor of 1.6. A direct comparison of our relaxation times to their results is problematic since they used a slightly different definition of the orientation correlation function. However, it is interesting to note that in the large-pore phase, the ratio τ lp /τ bulk = 2.4 found by Paesani 21 using classical MD agrees remarkably well with our value of 2.5.

Conclusions

The present work demonstrates that first-principles molecular dynamics simulations are a useful and computationally affordable tool for investigating structural and vibrational properties of metal-organic frameworks, providing information which is not easily accessible by experimental techniques. Our simulations revealed a markedly different hydration behavior of the narrow-and large-pore phases of MIL-53(Cr): In the narrowpore form of the material, water occupies proper crystalline sites, although with large vibrations around the equilibrium position. The large-pore phase, on the other hand, exhibits a more disordered hydration in which the adsorbed water shows a structure similar to that of bulk water, although with slightly different H bond arrangements. Regarding the dynamics of water in MIL-53(Cr), the mean square displacement of water molecules bound to µ 2 -OH groups in the framework is considerably larger in the narrow-than in the large-pore structure. This can be understood in terms of the steric hindrance due to other water molecules in the strongly hydrated largepore phase, which is absent in the narrow-pore case. Similarly, water reorientation dynamics are distinctly slower in the large-than in the narrow-pore phase. Moreover, reorientation of water molecules in the large-pore form is almost three times slower than in bulk water. This indicates that water in the pores of MIL-53(Cr), although structurally similar to bulk water, displays very clearly the impact of the confinement on its dynamics.

Fig. 1

 1 Fig. 1 Left: Structure of the hydrated narrow-pore phase, averaged over the trajectory. For clarity, the H atoms of the water molecules are not shown. Right: Structure of the hydrated large-pore phase, averaged over the trajectory. For clarity, the water molecules are not shown. Color code: Cr-brown, O-red, C-blue, H-white.

Fig. 2

 2 Fig.2Distribution of the angle between the µ 2 -OH group of the inorganic chain and the CrOCr plane, in the hydrated large-pore (lp) and narrow-pore (np) phases of MIL-53(Cr). Dashed lines represent Gaussian fits.

Fig. 3

 3 Fig. 3 Upper panel: time evolution of the O µ 2 -OH -O water distance of three water molecules near an µ 2 -OH group (see text for details). Lower panel: distribution of O µ 2 -OH -O water in the hydrated np and lp phases.

Fig. 4

 4 Fig. 4 Spectrum of the µ 2 -OH stretching vibration in the narrow-pore (black) and the large-pore (red) phase, obtained from Fourier-transforming the velocity autocorrelation function of the µ 2 -OH vibrations.

Fig. 5

 5 Fig. 5 Radial distribtion functions for O water -H µ 2 -OH in the narrowpore (black) and large-pore (red) phase.

Fig. 6

 6 Fig. 6 Snapshots of typical hydrogen bond arrangements of the µ 2 -OH-bound water molecule in MIL-53(Cr) lp phase.

Fig. 7

 7 Fig. 7 Radial distribution functions of water inside the MIL-53(Cr) large-pore phase (black). Upper panel: pair distribution of water oxygen (O water ), lower panel: pair distribution of O water and H. Red dashed line: bulk water.

Fig. 8

 8 Fig. 8 Left: Angle distribution for O water triplets in MIL-53(Cr), large-pore phase (black line). The red vertical line defines the tetrahedral angle of 109.5 • . The dashed red line represents bulk water. Right: Angle distribution for O water -H-O water triplets in MIL-53(Cr), large-pore phase (black line).

Fig. 9

 9 Fig. 9 Reorientation dynamics of water in the large-pore (upper panel) and narrow-pore (lower pannel) phase, characterized by the total correlation function C(t) (black line). The bulk water reference is shown as a dashed line. The colored curves represent component-wise correlation functions (eq. 2).

Acknowledgments

The authors acknowledge funding from the Agence Nationale de la Recherche under the project "SOFT-CRYSTAB" (ANR-2010-BLAN-0822) and generous allocation of computing time by GENCI (grant i2013087069).


François-Xavier Coudert

Daniela Kohen

Molecular Insight into CO 2 "Trapdoor" Adsorption in Zeolite Na-RHO

Among other zeolites, sodium-substituted zeolite Na-RHO exhibits a phenomenon known as "trapdoor" adsorption, with extra-framework cations acting as gates to small windows in the zeolite structures. While carbon dioxide can diffuse through these gates, methane and other guests cannot, leading to strong potential for gas separation. This effect has been attributed in the literature to specific cation-guest interactions which would allow to CO2 to pull the trapdoor open. We investigate here the gating phenomenon using ab initio molecular dynamics simulations combined with free energy methods. Our findings invalidate this previously-proposed mechanism, showing that presence of CO2 does not significantly affect the motion of the gating cation. We put forward and demonstrate an alternative mechanism, showing that the thermal motion of the cations are of large amplitude and that CO2 is able to squeeze while the gate is open, while nonpolar guests such as methane cannot. This brings an image of the mechanism that is closer to swinging doors than trapdoors.

INTRODUCTION

Effective carbon capture processes require materials with a high ability to discriminate between carbon dioxide and other gases. [START_REF] Cheung | Zeolites and related sorbents with narrow pores for CO2 separation from flue gas[END_REF] Microporous adsorbents, such as zeolites and metal-organic framework materials, have received a lot of attention in a drive to reduce greenhouse gas emissions due to their potential applications in gas separation, capture, and storage. These materials are also being studied given their potential applications in catalysis and other gas separation and storage processes. [START_REF] Keskin | Efficient Methods for Screening of Metal Organic Framework Membranes for Gas Separations Using Atomically Detailed Models[END_REF][START_REF] Yazaydin | Screening of Metal-Organic Frameworks for Carbon Dioxide Capture from Flue Gas Using a Combined Experimental and Modeling Approach[END_REF] A particularly promising class of materials for carbon capture are those exhibiting cation gating, a phenomenon which allows carbon dioxide but not other sorbents to permeate them. This is the case of some cationic zeolites whose extra-framework cations block the entrances of narrow pores connecting large cages, allowing cations to act as "trapdoors". Only carbon dioxide is then able to diffuse through the materials, bypassing the cations and moving from cage to cage. This cation gating gives rise to very high adsorption selectivities for CO2 over other small adsorbates, which cannot diffuse through the material. Some members of the Zeolite RHO and CHA (chabazite) family exhibit this behavior. Lozinska and coworkers [START_REF] Palomino | Zeolite Rho: a highly selective adsorbent for CO2/CH4 separation induced by a structural phase modification[END_REF][START_REF] Lozinska | Understanding Carbon Dioxide Adsorption on Univalent Cation Forms of the Flexible Zeolite Rho at Conditions Relevant to Carbon Capture from Flue Gases[END_REF][START_REF] Lozinska | Cation Gating and Relocation during the Highly Selective "Trapdoor" Adsorption of CO2 on Univalent Cation Forms of Zeolite Rho[END_REF] demonstrated that zeolite RHO materials with a Si/Al framework ratio around 4 have both good CO2 uptake and high CO2 selectivity with respect to small molecules like nitrogen, methane, and ethane. Similarly, chabazite zeolites CsCHA and KCHA exhibit remarkably high selectivity of carbon dioxide over methane. [START_REF] De Baerdemaeker | Gas Separation Trapdoors in Zeolites[END_REF][START_REF] Shang | Discriminative Separation of Gases by a "Molecular Trapdoor" Mechanism in Chabazite Zeolites[END_REF][START_REF] Shang | Determination of Composition Range for "Molecular Trapdoor" Effect in Chabazite Zeolite[END_REF] However, the molecular mechanism by which carbon dioxide passes through narrow pores blocked to other gases by cations is not well understood at the microscopic scale. Studies based on crystallographic and DFT evidence highlight the role of carbon dioxide/cation interaction and the collective effects requiring the presence of several CO2 molecules for the trapdoor to open. [START_REF] De Baerdemaeker | Gas Separation Trapdoors in Zeolites[END_REF][START_REF] Shang | Discriminative Separation of Gases by a "Molecular Trapdoor" Mechanism in Chabazite Zeolites[END_REF][START_REF] Shang | Determination of Composition Range for "Molecular Trapdoor" Effect in Chabazite Zeolite[END_REF] However, these approaches cannot take into account the dynamical nature of the phenomenon nor the realistic CO2 loading in the unit cell. Molecular dynamics simulations, on the other hand, offer an ideal tool to investigate the mechanism at play.

In this work we used ab initio molecular dynamics (AIMD) to study at the atomistic level the phenomenon of cation gating in the fully exchanged zeolite Na-RHO (Si/Al = 3.80). This Na-RHO structure is the most studied of the RHO family of zeolites, and the most promising for practical applications due to its selective adsorption properties and low cost. The zeolitic RHO topology [START_REF] Robson | Synthesis and Crystal-Structure of Zeolite Rho -New Zeolite Related to Linde Type-A[END_REF] (Figure 1) is typical of small-pore zeolites with window dimensions close to the kinetic diameter of the relevant gases and cage sizes that facilitate interaction with the adsorbing molecules. It has a 3-dimensional channel system composed of one-size cavities (αcages). Each α-cage is connected to six other α-cages by double 8ring pores (D8R). This connectivity gives rise to two interpenetrated pore systems, which are not interconnected. Na-RHO zeolite is known to exhibit some slight flexibility: when loaded with 1 bar of carbon dioxide the 8-rings are distorted, the α-cages become tetrahedral rather than cubic and the zeolite expands approximately 1.5% (but maintains its symmetry in the I43m space group). [START_REF] Lozinska | Understanding Carbon Dioxide Adsorption on Univalent Cation Forms of the Flexible Zeolite Rho at Conditions Relevant to Carbon Capture from Flue Gases[END_REF][START_REF] Lozinska | Cation Gating and Relocation during the Highly Selective "Trapdoor" Adsorption of CO2 on Univalent Cation Forms of Zeolite Rho[END_REF] This can be visualized in Figure 1, as well as in 3D using the CIF files for both structures which we have included as Supporting Information. Lozinska et al. have determined using X-ray diffraction that Na + cations preferentially occupy S8R (single 8-ring) sites and S6R (single 6-ring) sites (see Figure 1a), both when loaded with CO2 or not. The cations that sit in the S8R are the cations that block diffusion within the materials. The Si/Al ratio and site occupancies are such that all narrow pores are likely to blocked. Previous work on carbon dioxide within this zeolite using classical molecular dynamics, based on an empirical force field and a rigid description of the zeolite framework, suggested that electrostatic interactions between CO2 and zeolite produce a favorable binding site beyond the gate. This would opportunistically allow CO2 to partly enter the channel as the cation gate undergoes thermal fluctuations -although no actual crossing of the window was observed. [START_REF] Bamberger | Atomistic Simulations of CO2 during "Trapdoor" Adsorption onto Na-Rho Zeolite[END_REF] Figure 2 shows a rare event observed during classical simulation as the carbon dioxide molecule enters a D8R blocked by a cation (in blue). Presumably the fact that other guest molecules lack such strong interaction with the zeolite, result in their incapability to squeeze by the gating cation. These observations hint at a mechanism that differs from the one suggested earlier (see figure 1c) in that the crucial interaction is that of the carbondioxide/zeolite and thus the CO2 does not have an active role on moving the blocking Na + out of the way and into a different site. To get a good understanding of the mechanism for cation gating at a microscopic level, we need to go beyond the limitations of classical simulations or time-indepent experimental results. Given the inherent difficulty of obtaining experimental mechanistic information at the microscopic scale, we turned to high accuracy computational condensed matter simulations. Computational advances now make this type of system tractable for ab initio molecular dynamics (AIMD), and several studies have highlighted the potential of AIMD to shed light into framework-guest and cation-guest interactions in cationic zeolites. [START_REF] Coudert | Dipole moment, hydrogen bonding and IR spectrum of confined water[END_REF][START_REF] Coudert | Water nanodroplets confined in zeolite pores[END_REF][START_REF] Mace | Free energy barriers for CO2 and N-2 in zeolite NaKA: an ab initio molecular dynamics approach[END_REF] Below we present our results using AIMD to further demonstrate carbon dioxide's unique "gate opening behavior" in which its preference for binding inside the "gate" allows CO 2 to "squeeze" by the thermally fluctuating gate-keeping cation.

RESULTS AND DISCUSSION

Ab Initio Molecular Dynamics

We first performed unconstrained (free) ab initio molecular dynamics simulations of the Na-RHO zeolite, both without any guest molecules and loaded with 15 molecules of CO2 -which corresponds to a pressure of 1 atm of CO2 (See Computational Methods section for details). Figure 3 gives the reader a qualitative idea of the region of space sampled by each atom in the system. Figures 3a and 3b correspond to the same representative trajectory simulating a Na-RHO zeolite without any adsorbate, showing the average positions of all atoms (balls and cylinders) as well as trajectories of individual species (thin lines: Na + in Fig. 3a; framework atoms in

S6R S8R

Fig. 3b). Figure 3c corresponds to another representative trajectory simulating the Na-RHO zeolite loaded with carbon dioxide. We can see while framework atoms vibrate around their equilibrium positions, the sodium cations show relatively large amplitude of motion (up to 1 Å depending on the particular site) around the S8R and S6R sites, with locations that correspond well to the experimental crystallographic sites presented by Lozinska et. al. [START_REF] Lozinska | Cation Gating and Relocation during the Highly Selective "Trapdoor" Adsorption of CO2 on Univalent Cation Forms of Zeolite Rho[END_REF] (see SI). Finally, we note that carbon dioxide molecules sit preferentially on the S8R but also explore the cage, also in agreement with experimental evidence. In order to provide a quantitative analysis of atomic motion, and to look for vibrational modes that could correspond to a concerted motion involving the blocking cation and carbon dioxide, we analyzed the AIMD trajectories using the effective normal modes method. [START_REF] Martinez | Extracting effective normal modes from equilibrium dynamics at finite temperature[END_REF] The details of this analysis are shown in SI, and show that there is no unique vibration mode responsible for the gating phenomenon. Instead, we find that there are about 120 modes that invoke significant cation motion. With only 10 cations in the system, these modes also have contributions of many other atoms and correspond to heavily delocalized motions (e.g., "window breathing") that are a generic feature of zeolite frameworks. [START_REF] Coudert | Water nanodroplets confined in zeolite pores[END_REF] Given that a CO2 molecule entering the narrow pore or a cation moving from one site to another are likely to be rare events, it not surprising that there were not seen in the unconstrained simulations presented here. The computational demands of AIMD are such that free simulations long enough to witness one are prohibitive. To capture the mechanism we thus turned to free energy calculations using constrained molecular dynamics and the blue moon ensemble method [START_REF] Sprik | Free energy from constrained molecular dynamics[END_REF][START_REF] Carter | Constrained Reaction Coordinate Dynamics for the Simulation of Rare Events[END_REF] (see the Computational Methods section for details).

The mechanism governing the passage of carbon dioxide -unlike other gases, such as methane -from cage to cage in Na-RHO despite the presence of the blocking cation can be uncovered by evaluating the height of the free energy barriers involved. To circumvent the fact that an unconstrained MD simulation is very unlikely to sample a point high on the barrier, we use constrained ab initio molecular dynamics to "force" the system to cross the barrier. Following the blue moon approach, averaging the constraint forces gives a quantitative indication of how stable that particular configuration is, and integrating these forces along the reaction coordinate gives the configuration free energy as a function of constraint coordinate.

Free Energy Profile: Sodium cation leaving the blocking S8R site

Previous literature on cation gating in Na-RHO suggests [START_REF] Lozinska | Understanding Carbon Dioxide Adsorption on Univalent Cation Forms of the Flexible Zeolite Rho at Conditions Relevant to Carbon Capture from Flue Gases[END_REF][START_REF] Lozinska | Cation Gating and Relocation during the Highly Selective "Trapdoor" Adsorption of CO2 on Univalent Cation Forms of Zeolite Rho[END_REF] that the interaction of carbon dioxide with a blocking Na + cation allows it to "pull" the cation away briefly from the S8R site and into the S6R site. It is then suggested that during this brief interval CO2 enters the narrow pore. This scenario can be investigated by evaluating the free energy barrier for a cation to move from the blocking S8R site into the cage and onto a S6R site, in the presence and in the absence of carbon dioxide. For the suggested scenario to be confirmed the barrier must be significantly lower in the presence of guest molecules.

We thus compute free energy profiles as a function of the Na + cation coordinate perpendicular to a S8R plane, which corresponds to the direction of the cation leaving the S8R. This reaction coordinated, which we denote by x(COM-Na + ), correspond to the projection of the Na + cation's position on the x axis -as depicted in Figure 4a -while the cation is allowed to move freely in the yz plane. For convenience, its zero is chosen as the center of mass (COM) of the zeolitic framework.

The free energy profiles of the Na + cation along this direction, both with and without CO2, are plotted in Fig. 4b. They were obtained by running a series of constrained MD simulations with different values of x(COM-Na + ). Apart from the constrained cation, all other atoms in the system were allowed to move freely. It shows a Figure 4b shows the free energy profiles obtained in the absence and in the presence of 1 atm carbon dioxide. It is clearly visible that that the free energy barrier is not significantly smaller when there is carbon dioxide present (7.0 kcal/mol vs 6.7 kcal/mol). Thus, the calculated free energy profiles indicate that carbon dioxide molecules do not facilitate the movement of the Na + cation out of the S8R site. This is in contradiction with the mechanism proposed in earlier work, based on zero-Kelvin DFT calculations of the energy of cation displacement. This highlights the need to account for thermal effects and entropy, as done in our molecular dynamics simulations and through the use of free energy methods.

The free energy profiles do show however that, in the presence of CO2, Na + cations are less likely to migrate to the partner S8R (the ring that makes a D8R with the ring where the cation sits) which we denote S8R*. Moreover, the site they occupy is displaced by 0.5 Å into the cage (from x = 6.1 Å to x = 5.6 Å), and the free energy well becomes more symmetric. Figure 3c points to the reason behind this behavior: D8R contain a cation sitting on a S8R and a carbon dioxide sitting on S8R*, with an optimal Na + -CO2 distance "pushing" the cation slightly due to their interaction. This is confirmed by visual inspection of the constrained trajectories, highlighting the formation of CO2/Na + pairs in D8R.

It should finally be noted that the validity of these arguments rests on the appropriateness of the coordinate being used for the free energy profiles. When carbon dioxide is present the "true" reaction coordinate is likely to depend also on the positions of the adsorbed molecules. However the data presented here indicates that the mechanism is not likely to be that of carbon dioxide molecules interacting with a cation in a way that facilitates the opening the trapdoor. We thus turn to explore another option for the mechanism at play.

Free energy profile: CO2 molecule entering the S8R

From the mechanism hint gathered from classical simulations presented in the introduction (and depicted in Figure 2) and the free energy profile just presented, it appears worthwhile to study how the free energy profile of a guest molecule entering a S8R from the cage is affected by the presence of a nearby cation. Moreover, we note in Fig. 3a that Na + has a relatively large amplitude of motion, and from Fig. 4b we see that 2 kcal/mol allow a cation to move from its S8R site into the cage by ~1 Å. Thus at 298 K about 3% of the time the trapdoor is a little ajar, simply by means of thermal motion. Of course, a similar fraction of the time the cation would be exploring the space further from the cage, moving towards the partner S8R*. In the presence of carbon dioxide this other ring is likely to be occupied by another carbon dioxide and thus it is not likely that such an "exploration" would result in an opening for a passing adsorbate.

We thus calculated the free energy profiles for a carbon dioxide molecule entering the narrow pore from the cage, in two different situations: first, when a cation is sitting in its S8R site, with a constrained coordinate of x(COM-Na + ) = 5.76 Å; secondly, when the cation is displaced, with a constraint at x(COM-Na + ) = 4.68 Å. For both situations, the free energy profiles for CO2 were calculated as a function of the x-axis projection of the COM of the carbon dioxide (with the origin still set at the COM of the zeolitic framework). The free energy profiles for CO2 are plotted in Figure 5a. Immediately it is clear that a modest displacement of the sodium cation leads to a significant decrease, by a factor of 2, of the free energy barrier for carbon dioxide to enter the ring (from 11.6 kcal/mol to 6.2 kcal/mol). This decrease and the low value of the resulting free energy barrier, are key to explaining the observed diffusion of CO2 and the gating phenomenon. The large amplitude of the thermal motions of the Na + out of its crystallographic position allows the carbon dioxide to diffuse through.

The profiles in this case have been extended beyond the limit of the unit cell to allow for the fact that there is no symmetry between the two partners SR8. On the cage side there is a cation siting on (or A (kcal/mol)

x COM-Na + (Å) no CO 2 1 atm CO 2
x COM-Na+ near) the ring, while this is not the case on the other side. The minima at ~9 Å correspond to the stable S8R site for carbon dioxide. The fact that carbon dioxide has a stable site on the S8R is crucial to the understanding of the process. It is only because this is the case that the carbon dioxide is able to squeeze by the cation while the latter is out of its equilibrium position. Fig. 5c shows a free energy profile for carbon dioxide obtained when no cation was constrained so all Na + were allowed to evolve freely and none was occupying the relevant narrow pore. The profile clearly shows that the S8R is a favorable site for a carbon dioxide.

Moreover, visual inspection and comparison of the geometry of the framework when the CO2 is outside the S8R or constrained inside reveals that deformation of the framework, in the vicinity of the S8R, is necessary for the guest molecule to pass. The S8R aperture, measured in the narrowest dimension as the center-to-center distance between the two closest oxygen atoms, changes from an average value of 5.27 Å (when the CO2 is not in the S8R) to 5.57 Å (when the CO2 is passing through the S8R).

Finally, we tested our explanation of the mechanism by comparing the situation of carbon dioxide with that of methane, which is experimentally known not to permeate the material. We thus calculated the free energy profile for a methane molecule entering the S8R when the gate was partially open, i.e. with a Na + cation constrained at x(COM-Na + ) = 4.68 Å. Figure 6b shows that the barrier is significantly higher for methane than for carbon dioxide, in agreement with the experimental results. Furthermore, inspection of Fig. 5c, which was obtained when all cations were allowed to move freely and there was no cation in the D8R of interest, reveals that even in the absence of a cation the barrier for methane diffusion is rather large, much larger than that of carbon dioxide. Both findings are explained by the relatively low affinity of methane for the S8R, irrespective of the presence of a cation: methane sits preferably in the cage rather than in the rings, which are favorable for CO2 only because of electrostatic interactions.

Finally, we note that the diameter of the narrow pore within Na-RHO obtained using geometric criteria is only 1.9 Å. This is smaller than the kinetic diameter of both CO2 (3.3 Å in the narrowest dimension) and CH4 (3.8 Å). [START_REF] Ismail | Gas Separation Membranes: Polymeric and Inorganic[END_REF] Clearly, just looking at the size of the molecules and the pores is not a reliable indicator of the dynamical behavior of the system. The exclusion of methane cannot be accounted for merely by its size, but of its lack of favorable electrostatic interactions with the host material in the single-eight ring.

The difference in size of the adsorbate molecules, however, raises the possibility that size or shape selectivity play a role. a) b) c)

Figure 5. Free Energy Profile for a guest molecule entering the single-8 ring (S8R). a) Results for CO2 obtained when a cation is constrained to be at x(COM-Na + ) = 5.76 Å (S8R) in orange and when the cation is constrained x(COM-Na + ) = 4.68 Å in green. b) Profile when the cation is constrained at x(COM-Na + ) = 4.68 Å and thus the cation is out of its equilibrium position, carbon dioxide results are shown in green while those of methane are in purple. c) Profiles obtained when all cations can move freely (and none of them was located in a blocking position). Grey shows the methane results while aqua show carbon dioxide results.

CONCLUSIONS

We demonstrated here how the interaction of carbon dioxide with the zeolite Na-RHO's single eight ring (S8R) allows it to enter the narrow pore from the cage as the cation thermally fluctuates in and out of its blocking site. In this microscopic picture of the mechanism for selective CO2 adsorption, the sodium cations in the rings act as swinging doors agitated by thermal motion, through which the carbon dioxide can squeeze while they are ajar. Methane, lacking the strong adsorbate/zeolite interaction and thus an adsorption site at the gate entrance, experiences a much larger barrier even when the cation is a bit out of the gate. This is the root of the colossal selectivity reported experimental in earlier work.

As there is no experimental data providing direct proof of the diffusion mechanism, we used computational molecular simulation techniques to observe this mechanism. We relied on first principles molecular dynamics simulations of the fully flexible material, combined with free energy methods to obtain free energy profiles as a function of several relevant reaction coordinates -involving of the diffusing guest adsorbate and the blocking cation. To avoid the unfeasibly high computational cost of 2D free energy methods, we calculated cross-sections of this 2D free energy surface by combin- x(COM-ads) (Å)

CO 2 CH 4
ing constrained dynamics and the blue moon ensemble. This allowed us to reveal the microscopic mechanism at play.

This mechanism is schematized in the TOC. To move from one cage to another, a carbon dioxide approaches a narrow pore on the side of a S8R occupied by a cation that has thermally fluctuated a bit out of its site. This carbon dioxide squeezes into the pore via an activated process and then moves onto the partnering S8R*. Afterwards, the carbon dioxide can overcame the small barrier that separates the S8R* site from the cage on the right, and thus diffuse through the "blocked" the D8R.

This systematic study on zeolite Na-RHO underscores the power of ab initio molecular dynamics combined with free energy methods to study activated events such as guest diffusion in narrow pores blocked by cations. It is entirely transferable to other materials, with different zeolitic frameworks, different cations, or different adsorbate molecules. We expect this ab initio methodology to open the way to a more accurate description of guest diffusion at the microscopic level. Moreover, the main conclusions reached here on the importance of cation thermal motion and zeolite/carbon dioxide interactions shed a new light on previous models of narrow pore diffusion in zeolites. The mechanism observed is not specific to Na + or the RHO framework, and we expect that it could be observed for other cationic zeolites with narrow windows such as CHA or LTA. However, the quantitative details of the cation's thermal fluctuation motions, zeolite window size, and cation-CO2 interactions will in each case depend on window size, cation size and cation polarizability.

It puts in perspective the importance of the cation-carbon dioxide interactions, on which many theoretical studies focus but which may not always be a crucial factor limiting diffusion. We have shown that the cation trapdoor is very often a little open, even in the absence of guest molecules. Carbon dioxide molecules are unlikely to facilitate this behavior but they are likely to use the opportunity to get into the almost nominally blocked narrow pore, unlike other adsorbates that lack the preference for that narrow location. In the future, we intend to look at previously proposed diffusion mechanisms of polar species in narrow-pore zeolites to reanalyze them in light of our new findings. We also note that the molecular dynamics simulations and free energy studies performed in this work deal with the {adsorbent, adsorbate} system at or near thermodynamic equilibrium -whereas in industrial fixed-bed gas separation processes, kinetics might play an additional role in the permeation of small molecules. Molecular simulations at a higher scale, or mesoscopic modelling, would be required to study the influence of kinetics on dynamic separation of gas mixtures.

COMPUTATIONAL METHODS

The ab initio molecular dynamics simulations were performed using the CP2K simulation package, 19 which uses density functional theory to calculate atomic forces "on the fly". Its QUICKSTEP module uses the Gaussian and Plane Wave (GPW) method [START_REF] Lippert | A hybrid Gaussian and plane wave density functional scheme[END_REF] or a Gaussian and Augmented Plane Wave (GAPW) method, [START_REF] Iannuzzi | Inner-shell spectroscopy by the Gaussian and augmented plane wave method[END_REF] where the electronic density is expanded on a basis of plane waves. In the work presented here the use of the GAPW scheme for the O and Na atoms and the GPW scheme for Si, C and H atoms resulted in well-converged forces using a 300 Ry cutoff. We use the Perdew- Burke-Ernzerhof (PBE) functional [START_REF] Zhang | Comment on "Generalized gradient approximation made simple[END_REF] and D3 dispersion correction from Grimme. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] The core electrons are treated with the Goedecker-Teter-Hutter pseudopotentials. [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF] To avoid unphysical charges the sodium atoms were represented by DZVP MOLOPT SR [START_REF] Vandevondele | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases[END_REF] basis sets while all other atoms were represented by TZV2P [START_REF] Vandevondele | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases[END_REF] basis set. The nine outer electrons of Na were treated explicitly, but for all other atoms only their valence electrons were.

The motion of nuclei follows Newton's equation of motion with a time step of 0.5 fs. At each time step, the wave-function is optimized with the self-consistent field convergence criterion set to 1.0 × 10 -6 a.u. Simulations were performed in the NVT ensemble with a CSVR thermostat, [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] with a time constant of 100 fs. All simulations correspond to 298 K, the temperature at which carbon dioxide behavior within this material was studied experimentally.

Seven trajectories of at least 25 ps where obtained for an empty Na-RHO and seven others with 15 carbon dioxide molecules (the loading at 1 atm, the pressure at which experimental data is available ). In all cases one unit cell was simulated and periodic boundary conditions used. The initial positions for the cations and adsorbates used to start AIMD simulations were obtained using classical Grand Canonical Monte Carlo simulations while the framework atoms were those determined experimentally. The empty zeolite initial positions and those for the zeolite loaded with guest molecules are those of ref. [START_REF] Lozinska | Cation Gating and Relocation during the Highly Selective "Trapdoor" Adsorption of CO2 on Univalent Cation Forms of Zeolite Rho[END_REF] . The methods used are described in Bamberger [START_REF] Bamberger | Atomistic Simulations of CO2 during "Trapdoor" Adsorption onto Na-Rho Zeolite[END_REF] but briefly note that the interactions between carbon dioxide and Na-exchanged zeolite were modeled using the CCFF potential [START_REF] Fang | First principles derived, transferable force fields for CO2 adsorption in Na-exchanged cationic zeolites[END_REF] recently developed by Sholl and coworkers. This potential necessitates a rigid framework but allows the cations to move. More details on the initial configurations can be found on the SI.

As mentioned earlier, free energy profiles were obtained using constrained dynamics and the blue moon ensemble method. [START_REF] Sprik | Free energy from constrained molecular dynamics[END_REF][START_REF] Carter | Constrained Reaction Coordinate Dynamics for the Simulation of Rare Events[END_REF] To compute a free energy profile, at each point along the reaction coordinate the ensemble average force due to the constraint along the constraint direction was evaluated. When the constraint is a simple distance constraint like the one used in our work the average force can be related to the change in free energy by Here f is the force due to the constraint along the constraint direction and ζ* is the value of the constraint (fixed during a given simulation). The free energy as a function of the reaction coordinate (the free energy profile) is then obtained by integrating the average force along the reaction coordinate.

The constrained ab initio molecular dynamics simulations were performed using the CP2K code and the parameters detailed above, with constant temperature and fixed unit cell parameters. Note that CP2K computes f at each time step, as it is the Lagrange multiplier in the SHAKE algorithm used to constrain the dynamics.

To avoid strong oscillations in f the shake error was required to be less than 10 -8 a.u. In all cases the forces were averaged after 0.5 ps of evolution to capture the system after it had equilibrated. Errors reported correspond to standard deviation errors and were calculated using the block average method with 10 blocks (errors shown in the force graphs in the SI). A quadratic interpolation scheme was used to integrate the forces and thus to obtain the free energy profiles shown. The center of mass of the guest molecules was used to compute ζ* to avoid fluctuations related to rotations of the molecules. Care was taken to locate the constraint cation or molecule within the unit cell (avoiding the need to use the periodic boundary conditions in the evaluation of f). The initial configurations were chosen so the cations and guest molecules that were not constrained were relatively well equilibrated with the constrained object. More details are given in the SI.

ASSOCIATED CONTENT

Figure 1 .

 1 Figure 1. (a) The structure of zeolite Na-RHO when exposed to 1 bar of CO2. The framework is shown in blue. Note the 3-dimensional channel system composed of cavities (α-cages). Each α-cage is connected to six others by double-8 rings (D8Rs). Shown in grey and pink spheres are cages belonging to the two interpenetrated but not interconnected pore systems. (b) Representation of the Na + cation sites. Spheres showing preferred cation sites: S8R in orange and S6R in green. The figure also shows the D8R traced in blue and the S6R traced in black. Each unit cell has 6 D8R sites (and thus 12 S8R) and 8 S6R. The structure shown is that of the empty zeolite; comparing (a) and (b) note the slight distortion that happens upon carbon dioxide adsorption. (c) This figure shows the mechanism proposed by Lozinska and coworkers 5 at low carbon dioxide coverage. It has been reproduced with permission from their article. Note how in this mechanism the blocking cation first moves from one S8R site to its partner S8R site and then to a S6R site to allow for the passage of CO2. The green lines are meant to emphasize the role of the cation-carbon dioxide interaction (which Lozinska et al. presume would have an even larger stabilizing role in the presence of more carbon dioxide molecules).

Figure 2 .

 2 Figure 2. A hint of the microscopic mechanism for gating, observed from classical MD simulations: a carbon dioxide molecule entering a double 8-ring (D8R) blocked by a cation (in blue). Only the relevant adsorbate molecule and cation are shown in the figure (there are 9 more cations and 14 more guest molecules in the classical md simulation). The blocked D8R is highlighted.

Figure 3 .

 3 Figure 3. Representative ab initio molecular dynamics trajectories. a) System without guest molecules, with blue traces corresponding to the cation motion. For clarity only cations adsorbed on the central cage shown are displayed, and the rest of the zeolite is shown at its average positions. b) System without guest molecules, with traces corresponding to the motion of framework atoms; the cations are shown at their average positions. c) Na-RHO zeolite under 1 atm of carbon dioxide. For clarity only molecules adsorbed on the central cage shown are displayed; the rest of the zeolite is shown at its average positions.

Na 10 Figure 4

 104 Figure 4. a) Scheme of the reaction coordinate corresponding to a cation leaving the single 8-ring (S8R). For convenience the zero was chosen as the center of mass of the framework. Orange corresponds to the trajectory of a cation constrained to x(COM-Na + ) = 5.76 Å and green to a cation constrained at x(COM-Na + ) = 3.6 Å. The trace corresponds to the constrained trajectory (moving only in the yz plane) in the absence of adsorbed molecules. b) Free energy profile for a cation leaving the S8R site. Blue line: profile in the absence of carbon dioxide; red line: P(CO2) = 1 atm.

We thank Fabien Trousselet, Rodolphe Vuilleumier, Anne Boutin, and Marie-Laure Bocquet for stimulating discussions. We acknowledge access to HPC platforms provided by a GENCI grant (x2016087069).

Supporting Information

Cation and carbon dioxide locations. Effective normal modes. Free energy profiles. Details on the generation of initial configurations for AIMD. The Supporting Information is available free of charge on the ACS Publications website.

AUTHOR INFORMATION Corresponding Author

dkohen@carleton.edu ACKNOWLEDGMENT