Aoife C Fogarty

François-Xavier Coudert

Anne Boutin

Damien Laage

Reorientational dynamics of water confined in zeolites

Keywords:

We present a detailed molecular dynamics study of water reorientation and hydrogen-bond dynamics in a strong confinement situation, within the narrow pores of an all-silica Linde type A (LTA) zeolite. Two water loadings of the zeolite are compared with the bulk case. Water dynamics is found to be retarded in this extreme hydrophobic confinement and the slowdown is more pronounced at higher water loading. We show that water reorientation proceeds mainly by large-amplitude angular jumps, whose mechanism is similar to that determined in the bulk. The slowdown upon hydrophobic confinement is found to arise predominantly from an excluded-volume effect affecting the large fraction of water molecules lying at the interface with the zeolite matrix, with an additional minor contribution coming from a structuring effect induced by the confinement.

Introduction

Understanding the properties of water in confinement is of great importance in many contexts, including e.g. biochemistry where water can be found within narrow ion channels and protein cavities [1,2] , geology where water can be confined in mineral inclusions and in clays [3] , but also in silica nanopores [4][5][6] , Nafion membranes used for fuel cells [4,7] , carbon nanotubes [8] , metal-organic frameworks [9,10] and reverse micelles [11][12][13][14][15] .

Experiments and simulation studies have underlined the great impact of confinement on the structure, thermodynamics and dynamics of water [8,13,16] . The dynamics of confined water was measured to be retarded compared to the bulk situation, and several factors could be at the origin of this slowdown. However, their relative importance has not yet been clearly established. Possible factors include for example the chemical nature of the confining interface -hydrophilic or hydrophobic, polar or neutral [17] -the shape of the interface, the size of the cavity and the density and structure of the confined liquid.

Here, we focus on a model system to investigate the impact of extreme confinement. We consider water confined within an all-silica narrow-pore zeolite. Zeolites are microporous crystalline materials which have a very wide industrial use for hydrocarbon separation and in detergents [START_REF] Auerbach | Handbook of Zeolite Science and Technology[END_REF] . While most zeolites include sodium or aluminum ionic defects within their silica framework, we chose to study an all-silica zeolite [START_REF] Corma | [END_REF]20] . This zeolite contains cages of 11 Å in diameter which are hydrophobic due to the low polarity of the exposed oxygen atoms and to the absence of ionic or polar defects [21] , and in which water can be forced to enter under pressure [22] . Such a confining material is thus a good system to specifically investigate the effect of an extreme confinement on water dynamics, in the absence of specific interaction between water and the confining medium. We will focus here on water reorientation dynamics, which provides a good measure of the hydrogen(H)-bond network lability and dynamics, and which is accessible experimentally, for example via femtosecond infrared and NMR spectroscopies [23] .

The outline of the remainder of this paper is as follows. In Section 2, we present the key features of the zeolite material, the two water loadings under consideration and the simulation procedure. Section 3 describes the structure of the water nanopools within the zeolite cages, and Section 4 compares the water reorientation dynamics within the zeolite with the bulk reference. Section 5 then uses the extended jump model [23] recently developed by one of us to establish the molecular origin of the observed slowdown in water dynamics in strong hydrophobic confinement. Section 6 finally offers some concluding remarks and comparisons with other systems.

Methodology Zeolite Structure

The structure of the zeolite is the cubic Linde type A, or LTA. This structure contains α cages of diameter 13 Å, connected together in a cubic symmetry by 8-membered oxygen ring windows of diameter 6 Å, and smaller sodalite β cages of diameter 7 Å(see Figure 1). This all-silica LTA zeolite was synthesised and characterised experimentally [START_REF] Corma | [END_REF]20] . In our simulations, we employ a structure obtained by energy minimization at the Density Functional Theory level [24] . It is very similar to that determined in the experiments, containing a Si 24 O 48 unit cell with a cubic P m 3m symmetry and a lattice parameter of 11.9 Å. This unit cell contains one α supercage and one sodalite β cage.

The pore volume of the α cages was calculated with an approach similar to that developed by Connolly [25] and described in Refs. [26,27] . The solvent-accessible pore volume is delimited by the surface defined by rolling a probe sphere on the zeolite framework, defined as an assembly of hard spheres centered on the zeolite oxygens. The radii employed are 1.577 Å for the water oxygen probe and 1.5 Å for the zeolite oxygen sites, which correspond to half of the Lennard-Jones diameters of these respective sites within the force field employed here (see details below). The resulting accessible volume of one α supercage is V α =754 Å3 , which would yield a 5.6 Å radius for an ideal sphere. This value of the pore volume can then be used to determine the effective density of the confined liquid water. Two water loadings were simulated. In the first case (noted N=20 or high density), 20 water molecules per LTA unit cell were introduced, corresponding approximately to the saturation uptake in the liquid phase, as calculated by Grand Canonical Monte Carlo simulations [22] , and resulting in a 0.8 g.cm -3 mass density. This effective density is smaller than the bulk density in the same pressure and temperature conditions, consistent with earlier work on water density in other hydrophobic zeolites [26][27][28][29] . In the second case (N=15 or low density), 15 water molecules per unit cell were considered, corresponding to a 0.6 g.cm -3 mass density.

Simulations

Our simulations employ the rigid, non-polarisable SPC/E water model [30] which has been shown to properly reproduce the dynamics of bulk water at room temperature. [31] Due to the small size of the LTA unit cell, simulations are run on a 2 × 2 × 2 supercell to reduce the influence of periodic boundary conditions. The positions of the zeolite framework atoms are fixed. Water-zeolite interactions are modelled using the force field described in Ref. [22] , which was shown to correctly reproduce the adsorption thermodynamics of water in all-silica zeolites [28] . Long-range electrostatic interactions are treated by Ewald summation and both repulsion-dispersion and electrostatic interaction energies between the rigid zeolite framework and water molecules are precomputed on a grid (with a grid mesh of 0.2 Å) and stored for use during the simulation.

After equilibration in the canonical ensemble at 300 K, classical molecular dynamics are performed in the microcanonical ensemble with an average temperature of 304 K (±1 K). The propagation time step is 0.5 fs. All trajectories are run for more than 500 ps. The bulk water reference simulation with the same SPC/E water force field at 300 K is taken from a previous work [31] .

Jump analysis

The average jump mechanism, the jump angular amplitude and the jump time are determined using the same approach as described in detail in Ref. [31] .

Structure of water

We first characterize the structural properties of water confined inside the LTA zeolite. At the lower water loading (N=15), all water molecules lie within the α cages or within the windows connecting them. At the higher water loading (N=20), one water molecule per unit cell is found within the sodalite β cages while the remaining 19 water molecules are within the α cells. We therefore focus on water molecules within the α supercages.

We compute the radial distribution function between water oxygen sites. This is defined as

g(r) = 1 ρ n(r) V (r) , (1)
where ... is an average over all water oxygen sites, n(r) is the number of water oxygen atoms in a shell of width ∆r and radius r around the oxygen center under consideration, V (r) is the volume of that shell V (r) = 4πr 2 ∆r and ρ is the effective water oxygen number density within the α cage at the water loading under consideration, i.e. ρ = N/V α where V α =754 Å3 . The resulting radial distribution functions are presented in Fig. 2a. The comparison with the radial distribution function in the bulk shows that the location of the first peak, corresponding to the H-bond length, has not changed significantly upon confinement. In contrast, the second peak, which reveals the position of the second shell, seems to shift closer to the central water oxygen upon confinement.

Compared with the bulk situation, the height of the first peak appears to be enhanced in confinement at high density, and unchanged in confinement at low density. However, the radial distribution function was computed with the standard definition used for bulk liquids. In the zeolite cages, a large fraction of the cell volume is not accessible to water, and this fraction changes with distance. We therefore calculate a new radial distribution function better suited for a confined fluid,

g (r) = 1 ρ n(r) V (r) . (2)
Here V (r) is the volume of the shell at distance r from the oxygen site which is accessible to water. This pore volume is calculated with the Connolly approach described in Sec. 2, using a stochastic sampling of the shell V (r) to measure its overlap with the zeolite atoms and thus determine V (r).

The resulting renormalized radial distribution functions are presented in Fig. 2b. The accessiblevolume correction shows that the first peak is only slightly enhanced in confinement compared to the bulk situation. The peak is more enhanced at low water loading, possibly because each water molecule forms fewer H-bonds than at higher loading and than in the bulk, thus reinforcing the existing H-bonds. The second peak location is now similar in confinement and in the bulk, which shows that the apparent shift seen in the standard radial distribution function Fig. 2a was an artefact due to the variation of the accessible volume with the distance. The third peak is dramatically enhanced in confinement, especially at higher water loading. This peak occurs at a distance which is approximately the separation between water molecules which are at the zeolite interface and diametrically opposed within the same α cage. At this distance, the probability to find a pair of water molecules remains finite while the water-accessible volume V (r) becomes very small, leading to the observed enhancement of the radial distribution function.

Water reorientational dynamics

We now turn to the study of water reorientation and H-bond dynamics within the zeolite pores.

To measure the reorientation dynamics of water molecules, we compute the second-order reorientation time-autocorrelation function for a water OH-bond vector u, given by

C 2 (t) = P 2 [u(0) • u(t)] , (3)
where P 2 is the second-order Legendre polynomial. For each system, we calculate C 2 (t) as a function of each OH group's H-bond acceptor at the time origin. This H-bond acceptor can be either another water oxygen (O W ) or in a small fraction of cases a zeolite oxygen (O Z ). The following geometric H-bond criteria are used: R DA < 3.5 Å, R AH < 2.5 Å, and θ HDA < 30 • , where D is the water oxygen H-bond donor atom, A the H-bond acceptor oxygen atom and H a hydrogen atom. In these zeolite systems, the overwhelming majority of the water OH groups lie within the first hydration shell of the zeolite, defined as containing all OH groups whose oxygen lies within 5 Å of a zeolite oxygen. In both the high and low density systems, our results show that on average 96 % of the water OH groups are in the first hydration shell, of which on average 4 % point toward a zeolite oxygen atom and satisfy the H-bond geometric criteria. The reorientation time-correlation functions are shown in Figure 3a, along with that of bulk water for comparison. In all cases, C 2 (t) displays an initial rapid sub-ps decline followed by a longertime component. We now discuss in detail what an examination of these two time components can tell us about water dynamics in confinement in a zeolite system. All values referred to in the following discussion are summarised in Table 1.

Short time reorientation

We first focus on the initial rapid (< 200 fs) decay in C 2 (t). This sub-picosecond decay has been shown to be due to librational motions of the OH bond. [32] In an intact H-bond, the OH group performs a librational or wobbling motion in a cone whose axis is the vector between the water oxygen O W and the H-bond acceptor oxygen (either a zeolite or water oxygen) O A , within the frame defined by the O W -O A axis. This librational motion is the source of the initial rapid loss of correlation in the C 2 (t) function, and the size of the librational cone is a measure of H-bond strength, with cone angle increasing with decreasing H-bond strength. [32] We have calculated the cone angles for all possible H-bonds in our systems, following the procedure detailed in Ref. [32] .

We first consider the OH groups initially H-bonded to water oxygens O W . For these OH groups, the initial librational decay is approximately 20 % of the initial correlation (Fig. 3a), which corresponds to a cone semiangle of 25-27 • , a value very slightly below what is found for bulk water (27 • ). This suggests that the water-water H-bond is very slightly strengthened upon confinement. This is consistent with the small redshift of the OH stretch vibrational spectrum upon confinement that was computed for this zeolite [24] . The increase in water density, similarly, does not seem to have a strong impact on the cone semiangle.

Regarding the small fraction of water OH groups which point towards the zeolite and which engage in H-bonds with the zeolite oxygen atoms O Z , the initial librational decay leads to a ∼ 50 % loss of the orientational correlation (see Figure 3a), leading to cone semiangles of respectively 40 • and 35 • in the low and high water loadings. These values, which are markedly greater than for the O W H• • •O W H-bonds, show that the H-bond with the zeolite matrix is much weaker than with water. This is consistent with the large-amplitude librations already found for water next to sites whose polarity is low [33,34] . The polarity of the zeolite oxygen sites is so low that the water OH groups pointing toward the interface can almost be characterized as dangling OH bonds. These very weakly H-bonded groups appear as a small blue-shifted peak in the vibrational stretch spectrum [24,33] and have been extensively characterized experimentally at water-organic solvent interfaces [35] .

Long time reorientation

We now turn to the long-time component in the decay of the reorientation time-correlation function C 2 (t) (Eq. 3). It is quantified for each initial H-bond acceptor in each system via the secondorder reorientation time τ M D reor extracted from a monoexponential fit of C 2 (t) in the 2-10 ps time interval. The most striking observation is the slowdown in water reorientational dynamics caused by confinement, with reorientation times of 7.3 ps and 5.2 ps for water initially H-bonded to O W in the high and low density zeolite systems respectively, compared to 2.5 ps in bulk water. Secondly, increasing water density in the zeolite slows down reorientational relaxation, with reorientation times of 7.3 ps for an initial H-bond to O W and 5.5 ps for an initial H-bond to O Z at high density, compared to 5.2 ps and 4.8 ps respectively at low density. Finally, the identity of the initial Hbond acceptor also influences reorientational dynamics. At a given density, a water OH group initially H-bonded to O Z reorients more quickly than an OH group initially H-bonded to O W . The reorientation times are 5.5 ps vs 7.3 ps respectively in the high density system and 4.8 ps vs 5.2 ps respectively in the low density system.

In the next section, these observations will be explained using a theoretical model for water reorientation, the extended jump model, introduced by one of us. [23,36,37] This model will be used to precisely identify which molecular factors cause the measured slowdown in water dynamics within zeolites.

Extended jump model analysis

It has been shown that the main contribution to water reorientation both in the bulk and next to a range of solutes [23,38] , as well as in confinement in silica nanopores, [6] consists of large-amplitude angular jumps in which the reorienting OH group changes H-bond acceptors. We first demonstrate that these jumps are present within the zeolite pores and recall the main features of the extended jump model which connects the jump features to the water reorientation time. We then use the model to establish the relative importance of the different factors affecting water dynamics in strong hydrophobic confinement.

Characterization of angular jumps in zeolites

Large amplitude angular jumps were recently suggested [39] to occur in chabazite, a narrow-pore zeolite containing ionic defects. Here, we provide a confirmation and detailed characterization of these jumps.

We perform an a posteriori analysis of our MD trajectories, identifying each H-bond exchange event, i.e. each time a water hydroxyl group O*H* that was initially H-bonded to a water or zeolite oxygen A i becomes H-bonded to a different water oxygen A f . We define θ as the angle between the O*H* bond and the bisector plane of A i O * A f , and φ as the

A i O * A f angle. We also define distances R O * A i , R O * A f and R A i A f .
All parameters are illustrated in Figure 4. For each H-bond exchange event, we calculate these geometric parameters as a function of time during the 1 ps directly preceding and following the event. We examine more than 5000 such events for an initial O W acceptor and more than 1000 events for an initial O Z acceptor. We define a common time origin for each trajectory as the moment when θ = 0. This allows us to calculate the average time evolution of each geometric parameter as a function of the type of initial H-bond acceptor in each system, presented in Figure 5.

We find a large-amplitude angular jump mechanism for reorientation in each case, as evidenced by the data presented in Figure 5. The details of the mechanism are qualitatively similar in each zeolite system and in bulk water. The key features are the departure of the initial H-bond acceptor, seen in the increase in R O * A i after the transition state, the approach of the final acceptor to form a new H-bond, seen in the decrease in R O * A f , and the abrupt change in θ when the OH group switches from one acceptor to the other. We now discuss the quantitative differences between the systems.

Firstly, comparing the bulk water trajectories to those in the zeolite systems, we see that in confinement the average value of R O * A f before the H-bond exchange process is between 3.5 and 4.0 Å versus ∼ 4.5 Å in the bulk. Similarly, the final value of R O * Ai after the exchange is smaller in confinement than in the bulk. The first hydration shell around a given water molecule peaks at 2.8 Å in terms of oxygen-oxygen distance, while the second hydration shell lies between 4.4 and 4.6 Å (see Figure 2b). While in the bulk the new H-bond partner A f comes on average from the second hydration shell [31] , in confinement the new partner in a significant proportion of jumps comes from the first hydration shell, leading to a lower average initial value of R O * A f . A similar argument holds for the destination of the initial H-bond partner A i and the lower final value of R O * Ai in confinement.

Secondly, we see that R A i A f and φ are much larger for O Z →O W jumps than O W →O W jumps, with values of 3.75 to 5.5 Å and 70 to 110 • versus 3 to 3.5 Å and ∼ 60 • respectively. This is because a water OH group that initially points toward the zeolite must undergo a larger amplitude jump to find a new water H-bond acceptor which usually comes from behind it. A similar effect had been described for water at an apolar interface [33] [34] . The jump amplitude is defined as the value of the angle φ between the initial O*H*• • •A i and final O*H*• • •A f H-bond axes in the jump transition-state geometry. Again, values are given in Table 1. When the water density within the pore increases, the jump amplitude for O Z →O W jumps decreases slightly because the new partner can be more frequently found within the interfacial layer. The amplitude of the O W →O W jump is little affected by confinement in the zeolite or by a change in water density.

In summary, while confinement causes quantitative changes in the details of the H-bond exchange process, particularly for O Z →O W jumps, the fundamental jump mechanism underlying water reorientation remains unchanged.

Extended Jump Model

Having demonstrated the existence of large-amplitude angular jumps in the water reorientation mechanism in zeolites, we then go on to apply the extended jump model (EJM) [23] to connect these jumps to the measured water reorientation time in the zeolite systems.

Jumps are an activated process, and can usefully be treated as chemical reactions, i.e. a process whose rate constant can be determined by the free energy difference between its reactant and transition states. The jump kinetics can be followed via the cross time-correlation function (tcf) between stable states [31] I and F (in which the OH group forms a stable H-bond with the initial and final acceptors respectively). This jump tcf is given by [23] C jump (t

) = n I (0)n F (t) , (4)
where n I,F = 1 if the OH group is in stable state I and F respectively, and n I,F = 0 otherwise. Absorbing boundary conditions are used in the product state, so that only the first jump performed by each OH group is considered. The jump time, τ jump , which is the inverse of the rate constant for the H-bond acceptor exchange, can then be found via a monoexponential fit of 1 -C jump (t) in the 0.5-10 ps time interval. The amplitude of the jump angle ∆θ has an average bulk water value of 68 • , [31] and deviates from this value in the presence of a solute or in confinement.

A secondary and minor contribution to water reorientation arises from the tumbling of the molecular frame, i.e. the reorientation of an intact H-bond between jumps. [31] The time constant for this process, τ f rame reor , can be extracted from an exponential fit of the C 2 (t) reorientation correlation function of the O W -H vector between jumps, an approximation of the O W -O A vector between jumps.

The jump and frame-tumbling contributions to reorientation can then be combined in the extended jump model to give a predicted second-order reorientation time, thus [31] 1 τ EJM reor = 1

τ jump reor + 1 τ f rame reor , (5)
where the jump contribution is given by [31] τ jump reor = τ jump /[1 -sin(5∆θ/2)/(5 sin(∆θ/2)] .

The EJM model has been shown to give results in agreement with both NMR and pump-probe infrared experiments on bulk water, as well as with reorientation times calculated directly from molecular dynamics simulations. [37] It also holds true for water reorientation in confinement in zeolites, as we now show. The values resulting from application of the EJM are given in Table 1. Comparison of extended jump reorientation times from the model (τ EJM reor ) with reorientation times directly calculated from simulation (τ M D reor ) shows that the model correctly reproduces reorientation times in each case. The extended jump model can therefore be used to analyze and understand reorientational dynamics in these systems. We now use this model to elucidate the molecular origins of the effects of confinement and changes in density and the initial H-bond acceptor on reorientational dynamics.

Origin of slowdown

Our principal observation was the slowdown in water reorientation upon confinement in zeolites. We first concentrate on the comparison between reorientation times for OH groups initially Hbonded to water oxygen atoms, in either bulk water or the zeolite, since they represent the dominant population within the pores.

Table 1 shows that the frame contribution slows down dramatically in confinement. However, a comparison of the jump and frame reorientation times shows that the frame is a minor contribution to the overall reorientation time and it is in fact the slowdown in the jump component which is the main cause of the overall reorientation retardation. Within the EJM, the jump reorientation time is only determined by the jump amplitude and by the jump time, and the values reported in Table 1 unambiguously show that the τ reor retardation arises from a slowdown in the H-bond jump time τ jump .

We now analyze the origin of this slowdown in the jump time. The retardation factor compared with the bulk case is respectively 2.1 and 2.3 in the low and high zeolite water loading (see Table 1). A first cause of this slowdown comes from an excluded-volume effect for water OH groups lying at the interface with the zeolite wall. Our structural analysis has shown that they represent the overwhelming majority ( 96%) of the water molecules within the pores. Since part of the rate-limiting step in the H-bond jump mechanism is the approach of a new H-bond partner, the presence of a solute or an interface hinders this approach and leads to an excluded-volume effect at the jump transition state which results in a slowdown in the jump rate constant. This was described in detail for water molecules next to small hydrophobic convex solutes [40] and next to flat apolar interfaces [33,34] .

Starting from the simple anaytic prediction of the jump slowdown factor that was derived for water next to a spherical convex solute in Ref. [40] , it is straightforward to show that for water within a spherical concave pore, the slowdown factor can be approximated as

ρ jump = 1/π arccos r 2 -2rR + (R ‡ ) 2 -R ‡ d cos (∆θ) 2R ‡ sin (∆θ) (R -r) 2 -(d/2) 2 -1 , (7)
where R is the radius of the spherical cavity, r is the distance between the interfacial hydration layer and the minimum approach position next to the pore, R ‡ is the distance between oxygen atoms at the jump transition state, ∆θ is the amplitude of the angular jump and d is the average distance between two H-bonded water oxygen atoms. With typical values for these parameters (R ‡ = 3.2 Å, ∆θ = 60 • , d = 2.8 Å and r = 0.8 Å) and with the radius of the zeolite pores under consideration R = 5.6 Å, the predicted excluded-volume slowdown factor is ρ jump 1.9. This value is very close to what was previously determined next to a flat apolar surface [33,34] .

An additional contribution to the overall slowdown in the jump time arises from the enhancement of the water structure induced by the confining interfaces. While the effect of a single interface disappears after a few hydration layers, the confinement within the narrow pores of the zeolite studied here causes a strong structuring effect. This can be seen in the radial distribution functions presented in Fig. 2b which reveal an enhancement of the peaks and of the dips compared with the bulk situation. This increased local structure leads to a greater free energy cost for the elongation of the initial H-bond and for the approach of the new H-bond partner predominantly originating from the second shell, which were shown to bring a dominant contribution to the jump free energy barrier [31] . As shown in Ref. [31] , the jump free energy barrier can be approximated as

∆G ‡ w(R ‡ O * A i ) -w(R R O * A i ) + w(R ‡ O * A f ) -w(R R O * A f ) , (8)
where O * is the central rotating water oxygen, A i and A f are the initial and final H-bond acceptors, R ‡ is the distance at the jump transition state, and R R is the value in the reactants, i.e. the initial stable H-bond. w(r) is the potential of mean force along the oxygen-oxygen distance (Fig. 6), computed from the radial distribution function

w(r) = -k B T ln [g (r)] . (9)
Using the R ‡ and R R values determined from the average jump mechanisms in Fig. 5, the slowdown induced by the confinement-induced structuring is ρ = (∆G ‡ zeolite -∆G ‡ bulk )/k B T 1.2 ± 0.2. When combined together, these two slowdown factors lead to a satisfactory agreement with the directly computed retardation factor for the jump time (2.1 and 2.3 in the low and high density systems respectively). Our analysis has thus separated the contributions arising from an excludedvolume effect on interfacial water molecules and from a confinement-induced structuring effect. Although our estimates are not fully quantitative, our approach can already clearly show that the dominant effect arises from the steric excluded-volume effect due to the large fraction of water molecules lying at the interface, with an additional contribution from the structuring effect.

We further note that the faster reorientation dynamics measured at a given water density for the water OH groups which are initially pointing towards the zeolite oxygen atoms O Z compared to those H-bonded to water oxygens O W can also be rationalized by the extended jump model approach. Table 1 shows that this shorter reorientation time is due to a shorter jump time for H-bond exchange for O Z →O W jumps than for O W →O W jumps. This jump time in turn depends on the free energy cost to stretch the initial H-bond to the jump transition state [23] . As noted above from an examination of the librational cone size, an

O W H• • •O Z H-bond is weaker than an O W H• • •O W H-
bond, which facilitates the H-bond jumps and thus the reorientation.

Finally, we observe water dynamics to depend on the density of water in the zeolite. Increasing water density slows down reorientational relaxation. For an OH group initially H-bonded to a water oxygen, τ reor is 5.2 ps at low density and 7.3 ps at high density, while the corresponding values for an OH groups initially H-bonded to a zeolite oxygen are 4.8 and 5.5 ps. This slowdown is dominated by an increase in τ f rame reor , which is very sensitive to density. Frame tumbling is a diffusive motion, and is therefore retarded by the increasing viscosity which comes from increasing density. (We note that the frame tumbling time becomes infinite for OH groups H-bonded to O Z sites because the zeolite framework is fixed and this H-bond cannot tumble.)

Conclusions

We have shown that within the narrow, 11 Å-diameter hydrophobic pores of an all-silica zeolite, water reorientation occurs through large amplitude jumps. Water reorientation is retarded by a factor of 2-3 compared with the bulk situation. Most of this slowdown is due to an excluded-volume effect which affects the water molecules at the zeolite interface. Additional smaller contributions arise from the enhanced water structure induced by the confinement within the narrow zeolite pores and from the slower reorientation of intact H-bonds between successive H-bond jumps.

The present study has determined and analyzed the effect on water dynamics of a very strong confinement within a matrix which is hydrophobic and has no specific interaction with the water molecules. Other studies of wider nanopores [5,6] had found a retardation factor of approximately 2 for the interfacial water molecules, which is consistent with our present results. At the interface with large polar nanopores [5,6] and within reverse micelles of sizes similar to the zeolite cages studied here [11][12][13][14][15] , the presence of the strongly hydrophilic surfactant headgroups leads to a much greater slowdown of water dynamics. While the effect on water dynamics of hydrophilic defects within hydrophobic zeolites [21] will deserve a further specific study, these results together with our present study thus already suggest that the presence of the interface and its chemical nature have a much greater impact on water dynamics than confinement.

Aknowledgements

DL thanks Guillaume Stirnemann (ENS) for fruitful discussions.

Octogonal window

Sodalite cage

System

Low density water in zeolite High density water in zeolite Bulk water [31] 0.6 g.cm -3 0.8 g.cm 1: Characteristic times for water reorientation and the H-bond exchange process extracted directly from simulation and calculated via the extended jump model (EJM) [31] , as defined in the text. The librational cone angles are defined based on the reorientation time-correlation value after 200 fs using the procedure detailed in Ref. [32] .

Figure 7: Table of contents. Water dynamics in extreme confinement: the presence of the interface has a much greater impact on reorientational dynamics than the additional structuring of water due to confinement.

Figure 1 :

 1 Figure 1: Schematic view of the LTA structure showing an α cage (centre) with its eight surrounding sodalite β cages. α cages are connected together by octogonal windows.

Figure 2 :

 2 Figure 2: a) Standard radial distribution function (Eq. 1) between water oxygen atoms in bulk water and in the two zeolite systems. b) Same radial distribution functions corrected for accessible volume effects (Eq. 2).

Figure 3 :

 3 Figure 3: (a) Reorientation time-correlation functions C 2 (t) (Eq. 3) for water OH groups, in bulk water and in low and high density zeolite systems, as a function of the initial H-bond acceptor, which can be either another water oxygen (O W ) or a zeolite oxygen (O Z ). (b) Idem focussing on the initial decay. (c) Frame reorientation time-correlation functions, i.e. reorientation timecorrelation functions C 2 (t) for OH bonds between successive H-bond acceptor exchanges, for the same five subsets of OH groups.

Figure 4 :Figure 5 :

 45 Figure 4: Definition of the different geometric parameters used to follow the H-bond exchange event, as described in the text.

Figure 6 :

 6 Figure6: Potential of mean force as a function of water oxygen-water oxygen distance calculated using the g (r) radial distribution functions modified to take into account the volume accessible to the confined water (Eq. 2).

  -3 Jump type O W →O W O Z →O W O W →O W O Z →O W O W →O W

	τ MD reor (ps) τ frame (ps) reor Libration cone angle ( • )	5.2±0.1 14.2±0.1 27±1	4.8±0.1 ∞ 40±1	7.3±0.1 29.1±0.2 25±1	5.5±0.1 ∞ 35±1	2.5±0.1 5.4±0.1 27±1
	Jump angle ∆θ ( • )	58±2	103±2	60±2	74±2	68±2
	Jump time τ jump (ps)	6.9±0.2	5.5±0.2	7.6±0.3	5.5±0.2	3.3±0.1
	Jump reor. time τ jump reor (ps) EJM reor. time τ EJM reor (ps)	9.0 5.5±0.4	4.4 4.4±0.2	9.5 7.2±0.5	5.5 5.5±0.3	3.6 2.2±0.1
	Table


Giada Basile
email: basile@mat.uniroma1.it

AND Tomasz Komorowski
email: tkomorowski@impan.pl

Stefano Olla
email: olla@ceremade.dauphine.fr

Ceremade

DIFFUSION LIMIT FOR A KINETIC EQUATION WITH A THERMOSTATTED INTERFACE

Keywords: 2010 Mathematics Subject Classification. Primary, 60F17; Secondary, 82C70 Diffusion limit, thermostatted boundary conditions

We consider a linear phonon Boltzmann equation with a reflecting/transmitting/absorbing interface. This equation appears as the Boltzmann-Grad limit for the energy density function of a harmonic chain of oscillators with inter-particle stochastic scattering in the presence of a heat bath at temperature T in contact with one oscillator at the origin. We prove that under the diffusive scaling the solutions of the phonon equation tend to the solution ρ(t, y) of a heat equation with the boundary condition ρ(t, 0) ≡ T .

1. Introduction. We consider a linear phonon Boltzmann equation in contact with a heat bath at the origin. This equation describes the evolution, after a proper kinetic limit, of the phonon energy in a chain of harmonic oscillators with random scattering of velocities, where one oscillator is in contact with a heat bath at temperature T .

We denote by W (t, y, k) the energy density of phonons at time t ≥ 0, with respect to their position y ∈ R and frequency variable k ∈ T -the one dimensional circle, understood as the interval [-1/2, 1/2] with identified endpoints. The heat bath creates an interface localized at y = 0. Outside the interface the density satisfies ∂ t W (t, y, k) + ω′ (k)∂ y W (t, y, k) = γ 0 LW (t, y, k), (t, y, k

) ∈ R + × R * × T * , (1.1)
Here R + := (0, +∞), R * := R \ {0}, T * := T \ {0},

The parameter γ 0 > 0 represents the phonon scattering rate. The scattering operator L, acting only on the k-variable, is given by

LF (k) := T R(k, k ′ ) [F (k ′ ) -F (k)] dk ′ , k ∈ T
for F belongs to B b (T) -the set of bounded measurable, real valued functions. Here

ω′ (k) = ω ′ (k) 2π , k ∈ T, (1.2)
where ω : T → [0, +∞) is the dispersion relation of the harmonic chain.

The interface conditions that describe the interaction of a phonon with a thermostat (placed at y = 0), at temperature T > 0, are given as follows:

-the outgoing densities are given in terms of the incoming ones as W (t, 0 + , k) = p -(k)W (t, 0 + , -k) + p + (k)W (t, 0 -, k) + T g(k), for 0 < k ≤ 1/2, W (t, 0 -, k) = p -(k)W (t, 0 -, -k) + p + (k)W (t, 0 + , k) + T g(k), for -1/2 < k < 0.

( 1.3) where p -, p + , g : T → (0, 1] are continuous and

p + (k) + p -(k) + g(k) = 1. (1.4)
In other words, p -(k) and p + (k) are the reflection and transmission coefficients across the interface, respectively. They correspond to the probabilities of the phonon being reflected, or transmitted, by the interface.

The quantity T g(k) is the phonon production rate by the thermostat as well as the absorption rate of the frequency k phonon by the interface. The parameter T > 0 is the heat bath temperature. This equation has been obtained in [START_REF] Basile | Wigner functions and stochastically perturbed lattice dynamics[END_REF], without the heat bath, as the Boltzmann-Grad limit of the energy density function for a microscopic model of a heat conductor consisting of a one dimensional chain of harmonic oscillators, with inter-particle scattering conserving the energy and volume. In the presence of the thermostat, but with no scattering (the case γ 0 = 0), the limit has been proved [START_REF] Komorowski | High frequency limit for a chain of harmonic oscillators with a point Langevin thermostat[END_REF]. It is believed that the limit also holds in case of the presence of scattering, i.e. when γ 0 > 0.

We are interested in the asymptotics of the solutions to (1.1) under the diffusive scaling, i.e. the limit, as ǫ → 0, for W ε (t, y, k) = W (t/ε 2 , y/ε, k), i.e. we consider the equation

∂ t W ε (t, y, k) + 1 ε ω′ (k)∂ y W ε (t, y, k) = γ ε 2 T R(k, k ′ ) [W ε (t, y, k ′ ) -W ε (t, y, k)] dk ′ , y = 0, W ε (0, y, k) = W 0 (y, k), (1.5)
with the interface conditions (1.3). Let R(k) = R(k, k ′ )dk ′ . In our main result, see Theorem 2.2 below, we prove that under the assumption

T ω ′ (k) 2 R(k) dk < +∞, (1.6)
and some other technical hypotheses, formulated in Sections 2.2 and 2.3 below, for any

G ∈ C ∞ 0 (R × T) -compactly supported, smooth function -we have lim ǫ→0 R×T W ε (t, y, k)G(y, k)dydk = R×T ρ(t, y)G(y, k)dydk, (1.7)
where

∂ t ρ(t, y) = D∂ 2 y ρ(t, y), (t, y) ∈ R + × R * , (1.8)
ρ(t, 0 ± ) ≡ T (1.9)

ρ(0, y) = ρ 0 (y) := T W 0 (y, k)dk.
The diffusion coefficient is given by

D = 1 γ ω ′ (k)(-L) -1 ω ′ (k) dk (1.10)
that is finite by the assumption (1.6) and the properties of R(k, k ′ ) made in section 2.2.

The result implies that only the absorption and the creation of phonons at the interface matter in this diffusive scale. Phonons that are reflected or transmitted will come back to the interface, due to scattering, and eventually get absorbed in a shorter time scale.

The diffusive limit has been considered, without the presence of interface, in [START_REF] Basile | From a kinetic equation to a diffusion under an anomalous scaling[END_REF][START_REF] Jara | Limit theorems for additive functionals of a Markov chain[END_REF][START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF]. It has been shown

there that, if (1.6) is in force, then the solutions of the initial problem (1.5) converge, as in (1.7), to ρ(t, y)

-the solution of the Cauchy problem for the heat equation (1.8). When the condition (1.6) is violated a superdiffusive scaling may be required and the limit could be a fractional diffusion. This case has been also considered in [START_REF] Basile | Convergence of a kinetic equation to a fractional diffusion equation, Markov Process[END_REF][START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations: a Hilbert expansion approach[END_REF][START_REF] Jara | Limit theorems for additive functionals of a Markov chain[END_REF][START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF].

The case of the diffusive limit of the solution of a kinetic equation with an absorbing boundary has been considered in e.g. [START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF][START_REF] Bensoussan | Boundary Layers and Homogenization of transport processes[END_REF][START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF][START_REF] Bardos | Diffusion approximation and computation of the critical size[END_REF][START_REF] Bardos | The Diffusion Approximation for the Linear Boltzmann Equation with Vanishing Scattering Coefficient[END_REF]. The diffusive limit with some other boundary conditions has also been discussed in the review paper [START_REF] Bardos | Half-Space Problems for the Boltzmann Equation: A Survey[END_REF], see the references contained therein. A related result for the radiative transport equation with some reflection/transmission condition has been obtained in [START_REF] Bal | Diffusion approximation of radiative transfer problems with interfaces[END_REF] for the steady state, giving rise to different boundary conditions. We are not aware of a similar result in the dynamical case, as considered in the present paper.

The result for a fractional diffusive limit with the interface condition is a subject of the paper [START_REF] Komorowski | Fractional diffusion limit for a kinetic equation with an interface -probabilistic approach[END_REF].

2. Preliminaries and the statement of the main result.

W (t, 0 + , k) = p -(k) W (t, 0 + , -k) + p + (k) W (t, 0 -, k)), for 0 < k ≤ 1/2, ( 2.1)
and

W (t, 0 -, k) = p -(k) W (t, 0 -, -k) + p + (k) W (t, 0 + , k), for -1/2 < k < 0. (2.2) Definition 2.1. A function W : R+ × R × T → R is called a (weak) solution to equation (1.1) with the interface conditions (2.1) and (2.2), provided that it belongs to L 2 loc (R + , L 2 (R × T)), its restrictions W ι to R + × R ι × T, ι ∈ {-, +} extend to continuous functions on R + × Rι × T that satisfy (2.1) and (2.
2), and

0 = +∞ 0 R×T W (t, y, k) [∂ t ϕ(t, y, k) + ω′ (k)∂ y ϕ(t, y, k) + γLϕ(t, y, k)] dtdydk + R×T W 0 (y, k)ϕ(0, y, k)dydk (2.3) for any test function ϕ ∈ C ∞ 0 ( R+ × R * × T), 2.
2. Assumption about the dispersion relation and the scattering kernel. We assume that ω(•) is even, belongs to C ∞ (T \ {0}), i.e. it smooth outside k = 0. Furthermore we assume that ω(•) is unimodal,

that implies that kω ′ (k) ≥ 0 for k ∈ (-1/2, 1/2).
We assume that the scattering kernel is symmetric

R(k, k ′ ) = R(k ′ , k), (2.4)
positive, except for 0 frequency, i.e.

R(k, k ′ ) > 0, (k, k ′ ) ∈ T 2 * (2.5)
and the total scattering kernel

R(k) := T R(k, k ′ )dk ′ (2.6)
is such that the stochastic kernel

p(k, k ′ ) := R(k, k ′ ) R(k) ∈ C ∞ (T 2 ).
(2.7)

In addition we assume that (1.6) is in force.

Example. Suppose that

R(k) ∼ R 0 | sin(πk)| β , |k| ≪ 1
for some β ≥ 0 and R 0 > 0 and

ω ′ (k) ∼ 2ω ′ 0 sign(k) | sin(πk)| κ , |k| ≪ 1 (2.8)
for some κ ≥ 0. Then (1.6) holds, provided that 0 ≤ β < 1 + 2κ.

(2.9) 2.3. About the reflection, transmission and absorption coefficients. In [START_REF] Komorowski | High frequency limit for a chain of harmonic oscillators with a point Langevin thermostat[END_REF] the coefficient p ± (k) and g(k) are obtained from the microscopic dynamics and depends on the dispersion relation as follows.

Let γ > 0 (the thermostat strength) and g(λ

) := 1 + γ T λdk λ 2 + ω 2 (k) -1
, Re λ > 0.

(2.10)

It turns out, see [START_REF] Komorowski | High frequency limit for a chain of harmonic oscillators with a point Langevin thermostat[END_REF], that

|g(λ)| ≤ 1, λ ∈ C + := [λ ∈ C : Re λ > 0]. (2.11)
The function g(•) is analytic on C + . By Fatou's theorem we know that

ν(k) := lim ǫ→0+ g(ǫ -iω(k)) (2.12)
exists a.e. in T and in any L p (T) sense for p ∈ [1, ∞). Denote

g(k) := γ|ν(k)| 2 |ω ′ (k)| 2 , P(k) := γν(k) 2|ω ′ (k)| p + (k) := |1 -P(k)| 2 (2.13) p -(k) := |P(k)| 2 .
It has been shown in Section 10 of [START_REF] Komorowski | High frequency limit for a chain of harmonic oscillators with a point Langevin thermostat[END_REF] that

Re ν(k) = 1 + γ 2|ω ′ (k)| |ν(k)| 2 . (2.14)
This identity implies in particular that (1.4) is in force.

2.4. Scaled kinetic equations and the formulation of the main result. Consider W ε the solution of a rescaled kinetic equations (1.5). Our main result can be stated as follows.

Theorem 2.2. Suppose that W 0 (y, k) = T + W 0 (y, k), where W 0 ∈ L 2 (R × T). Under the assumptions made

about the scattering kernel R(•, •) and dispersion relation ω(•), for any test function ϕ ∈ C ∞ 0 (R + × R × T) we have lim ǫ→0 +∞ 0 dt R×T W ε (t, y, k)ϕ(t, y, k)dydk = +∞ 0 dt R×T ρ(t, y)ϕ(t, y, k)dydk, (2.15)
where ρ(•, •) is the solution of the heat equation

∂ t ρ(t, y) = D∂ 2 y ρ(t, y), (t, y) ∈ R + × R * , (2.16
)

ρ(0, y) = ρ 0 (y) := T W 0 (y, k)dk, ρ(t, 0) ≡ T, t > 0.
Here, the coefficient D > 0 is given by (4.6) below.

Defining W ε = W ε -T , one can easily see that it also satisfies (1.5) with and

W ε (t, 0 + , k) = p -(k) W ε (t, 0 + , -k) + p + (k) W ε (t, 0 -, k), for 0 ≤ k ≤ 1/2, (2.17
W ε (t, 0 -, k) = p -(k) W ε (t, 0 -, -k) + p + (k) W ε (t, 0 + , k), for -1/2 ≤ k ≤ 0. ( 2

.18)

The initial condition W ε (0, y, k) = W 0 (y, k) := W 0 (y, k) -T belongs to L 2 (R × T). This means that it is enough to prove Theorem 2.2 for T = 0. This proof is presented in Section 4.

3. Some auxiliaries.

3.1. Some functional spaces. Let H 1 + be the Hilbert obtained as the completion of the Schwartz class S(R × T) in the norm

ϕ 2 H 1 + := ϕ 2 L 2 (R×T) + +∞ 0 T |ω ′ (k)|[∂ y ϕ(y, k)] 2 dydk
Similarly, we introduce H 1 -. Let H be the Hilbert space obtained as the completion of S(R × T) in the norm

ϕ 2 H := ϕ 2 H 1 - + ϕ 2 H 1 + . (3.1)
Let also

ϕ 2 H0 := T |ω ′ (k)|g(k) [ϕ(0 + , k)] 2 + [ϕ(0 -, k)] 2 dk + R×T 2 [ϕ(y, k ′ ) -ϕ(y, k)] 2 dydkdk ′ 3.2. Apriori bounds.
Computing the time derivative we have 1 2

d dt W ε (t) 2 L 2 = - γ 2ε 2 ∞ -∞ dyD( W ε (t, y, •)) - 1 2ε T ω′ (k) W ε (t, 0 -, k) 2 -W ε (t, 0 + , k) 2 dk, (3.2)
with

D(f ) := T 2 R(k, k ′ )[f (k) -f (k ′ )] 2 dkdk ′ . (3.3)
Taking into account (2.17) and (2.18) we obtain

T ω′ (k) [ W ε (t, 0 -, k)] 2 -[ W ε (t, 0 + , k)] 2 dk = 1/2 0 ω′ (k) [ W ε (t, 0 -, k)] 2 -p -(k) W ε (t, 0 + , -k) + p + (k) W ε (t, 0 -, k) 2 dk + 0 -1/2 ω′ (k) p -(k) W ε (t, 0 -, -k) + p + (k) W ε (t, 0 + , k) 2 -[ W ε (t, 0 + , k)] 2 dk.
After straightforward calculations (recall that coefficients p ± (k) are even, while ω′ (k) is odd) we conclude that the right hand side equals

1/2 0 ω′ (k) W ε (t, 0 -, k) 2 + W ε (t, 0 + , -k) 2 1 -p 2 + (k) -p 2 -(k) -4p -(k)p + (k) W ε (t, 0 + , -k) W ε (t, 0 -, k) dk. Since p + (k) + p -(k) ≤ 1 we have 1 -p 2 + (k) -p 2 -(k) ≥ 0. In addition, det      1 -p 2 + (k) -p 2 -(k) -2p -(k)p + (k) -2p -(k)p + (k) 1 -p 2 + (k) -p 2 -(k)      = 1 -(p + (k) + p -(k)) 2 1 -(p + (k) -p -(k)) 2 .
Using (1.4) we conclude that the quadratic form

(x, y) → 1 -p 2 + (k) -p 2 -(k) (x 2 + y 2 ) -4p -(k)p + (k)xy is positive definite as long as p + (k) + p -(k) < 1.
The eigenvalues of the form can be determined from the equation

0 = 1 -λ -(p + (k) + p -(k)) 2 1 -λ -(p + (k) -p -(k)) 2 ,
which yields

λ + := 1 -(p + (k) -p -(k)) 2 , λ -:= 1 -(p + (k) + p -(k)) 2 and λ + > λ -. Note that 2g(k) ≥ λ -= g(k) [1 + p + (k) + p -(k)] ≥ g(k).
Equality (3.2) allows us to obtain the following apriori bounds

W ε (t) 2 L 2 (R×T) ≤ W 0 2 L 2 (R×T) , t 0 ds R D( W ε (s, y, •))dy ≤ ε 2 γ W 0 2 L 2 (R×T) , t 0 ds 1/2 0 ω′ (k)g(k)dk W ε (s, 0 -, k) 2 + W ε (s, 0 + , -k) 2 ≤ ε W 0 2 L 2 (R×T) .
(3.4) By (2.17) and (2.18) we obtain that

W ε (s, 0 + , k) 2 ≤ W ε (s, 0 -, k) 2 + W ε (s, 0 + , -k) 2 , k ∈ (0, 1/2) W ε (s, 0 -, k) 2 ≤ W ε (s, 0 -, -k) 2 + W ε (s, 0 + , k) 2 , k ∈ (-1/2, 0). (3.5)
Then using the unimodality of ω(k) it follows that

t 0 ds T dk|ω ′ (k)|g(k) W ε (s, 0 -, k) 2 + W ε (s, 0 + , k) 2 ≤ 2ε W 0 2 L 2 (R×T) . (3.6)
3.3. Uniform continuity at y = 0. Suppose that y > 0. Let

V ε (t, y, k) := t 0 W ε (s, y, k)ds. (3.7) Since W ε (s, y, k) satisfies (1.5) we can write ε W ε (t, y, k) -W ε (0, y, k) + ω′ (k)∂ y V ε (t, y, k) = F ǫ (t, y, k), (3.8)
with

F ǫ (t, y, k) := γ ε T R(k, k ′ ) [V ε (t, y, k ′ ) -V ε (t, y, k)] dk ′ , y = 0.
Hence, using Cauchy-Schwarz inequality, we get

F ǫ (t, •) 2 L 2 (R×T) = γ ε 2 R×T dydk t 0 ds T R(k, k ′ ) W ε (s, y, k ′ ) -W ε (s, y, k) dk ′ 2 ≤ γ ε 2 R×T dydk t 0 ds T R(k, k ′ )dk ′ t 0 ds T R(k, k ′ ) W ε (s, y, k ′ ) -W ε (s, y, k) 2 dk ′ ≤ t R ∞ γ ε 2 t 0 ds R D( W ε (s, y, •))dy.
Using the second estimate of (3.4) we conclude that for each t 0 > 0 sup

ǫ∈(0,1] sup t∈[0,t0] F ǫ (t, •) L 2 (R×T) ≤ γt 0 R ∞ W 0 2 L 2 (R×T) < +∞. (3.9)
From (3.8) we conclude that

∂ y V ε (t, y, k) = Fǫ (t, y, k) ω ′ (k) , y, ω′ (k) = 0, (3.10)
where

Fǫ (t, y, k) = F ǫ (t, y, k) -ε W ε (t, y, k) -W ε (0, y, k) .
From (3.9) and the first estimate of (3.4) we conclude sup

ǫ∈(0,1] sup t∈[0,T ] Fǫ (t, •) L 2 (R×T) =: F * (T ) < +∞.
We have

V ε (t, •) H 1 ± ≤ Fǫ (t, •) L 2 (R×T) t ≥ 0. (3.11)
Since Vǫ (t) = W ǫ (t), from the first estimate of (3.4) we conclude that for any t 0 > 0 sup

ǫ∈(0,1] Vε L ∞ ([0,t0];L 2 (R×T))
< +∞.

From (3.11) we get also (cf (3.1))

sup ǫ∈(0,1] V ε L ∞ ([0,t0];H) < +∞.
Summarizing, we have shown the following.

Proposition 3.1. For any t 0 > 0

C(t 0 ) := sup ǫ∈(0,1] V ε L ∞ ([0,t0];H) + Vε L ∞ ([0,t0];L 2 (R×T)) < +∞ (3.12)
and

lim ε→0+ V ε L ∞ ([0,t0];H0) = 0. Denote by W 1,∞ 0 ([0, t 0 ]; L 2 (R × T)) the completion of the space of smooth functions f : [0, t 0 ] → L 2 (R × T)
satisfying f (0) = 0, with respect to the norm

f W 1,∞ 0 ([0,t0];L 2 (R×T)) := sup t∈[0,t0] ḟ L 2 (R×T) .
As a consequence of the above proposition we immediately conclude the following.

Corollary 3.2. The family

(V ε (•)) ǫ∈(0,1] is bounded in W 1,∞ 0 ([0, t 0 ]; L 2 (R× T))∩L ∞ ([0, t 0 ]; H) for aby t 0 > 0. Any ⋆-weak limit point V (•) of V ε (•), as ε → 0+, satisfies the following: 1) V (t, y, k) ≡ V (t, y) := T V (t, y, k)dk for (t, y, k) ∈ R + × R × T) 2) the mapping R + × R ι ∋ (t, y) → V (t, y) extends to a mapping from C( R+ × Rι ), ι ∈ {-, +}, 3) V (t, 0 ± ) = 0 for each t > 0, 4) V (0, y) ≡ 0, y ∈ R.
4. Proof of Theorem 2.2. Thanks to the above estimates, we conclude that the solutions W ε (•) are ⋆-

weakly sequentially compact in L ∞ [0, +∞); L 2 w (R × T) , where L 2 w (R × T) denotes L 2 (R × T)
with the weak topology. Suppose that W (t, y, k) is a limiting point for some subsequence ( W εn (s, y, k)), where ε n → 0. For convenience sake we shall denote the subsequence by ( W ε (s, y, k)). Thanks to (3.4) for each t > 0 we have

(cf (3.3)) lim ǫ→0 t 0 ds R D W ε (s, y, •) dy = 0 (4.1)
thus W (t, y, k) ≡ ρ(t, y), for a.e. (t, y, k

) ∈ R + × R × T. Lemma 4.1. Equation -LX 1 = ω′ (4.2)
has a unique solution such that

T X 1 (k)R(k)dk = 0 and T X 2 1 (k)R(k)dk < +∞. (4.3)
Proof. Let µ be a Borel probability measure on T given by

µ(dk) = R(k) R dk, where R := T R(k)dk.
We can reformulate (4.2) as

X 1 -P X 1 = ω′ R , (4.4)
where, by virtue of (1.6), the right hand side belongs to L 2 (µ) and P : L 2 (µ) → L 2 (µ) is a symmetric operator on L 2 (µ) given by

P F (k) := T p(k, k ′ )F (k ′ )dk ′ , F ∈ L 2 (µ).
The operator is a compact contraction and, since

T F (k)(I -P )F (k)F (k)µ(dk) = D(F )
we conclude that 1 is a simple eigenvalue, with the respective eigenspace spanned on the eigenvector F 0 ≡ 1.

Thus the conclusion of the lemma follows, as ω′ /R ⊥ F 0 . Proof. The claim made in (4.6) follows immediately from the fact that X 1 = 0 and D = D(χ).

To prove (4. where ϕ ∈ C ∞ 0 (R + × R * ) and χ j (t, y, k), j = 1, 2 are yet to be determined. We can write

0 = +∞ 0 R×T ∂ t W ε (t, y, k)ϕ ε (t, y, k) dtdydk = +∞ 0 R×T W ε (t, y, k)∂ t ϕ ε (t, y, k) + ∂ t W ε (t, y, k)ϕ ε (t, y, k) dtdydk = +∞ 0 R×T W ε (t, y, k)∂ t ϕ ε (t, y, k) - 1 ε ω′ (k)∂ y W ε (t, y, k)ϕ ε (t, y, k) + γ ε 2 L W ε (t, y, k)ϕ ε (t, y, k) dtdydk. (4.8)
Since the support in the y variable of the test function ϕ ǫ (t, y, k) will turn out to be isolated from 0, the integration by parts yield the equation

0 = +∞ 0 R×T W ε (t, y, k) ∂ t ϕ ε (t, y, k) + 1 ε ω′ (k)∂ y ϕ ε (t, y, k) + γ ε 2 Lϕ ε (t, y, k) dtdydk. (4.9)
Substituting from (4.7) we obtain that the term in the brackets has the form

1 ǫ I + II + ǫIII ǫ , (4.10)
with

I := γLχ 1 + ω′ (k)∂ y ϕ(t, y), II := γLχ 2 + ω′ (k)∂ y χ 1 + ∂ t ϕ(t, y), (4.11)
III ǫ := ∂ t χ 1 (t, y, k) + ε∂ t χ 2 (t, y, k) + ω′ (k)∂ y χ 2 .
We stipulate that I = 0 and II = ∂ t ϕ(t, y) + D∂ 2 yy ϕ(t, y). (4.12)

The first condition yields Since (D + ω′ X 1 /γ)/R belongs to L 2 (µ) and its orthogonal to constants, we can solve the equation

χ 1 (t, y, k) = - 1 γ ∂ y ϕ(t, y)L -1 ω′ (k) = - 1 γ ∂ y ϕ(t, y)X 1 (k), ( 4
LX 2 = D + 1 γ ω′ X 1
using the same argument as in Lemma 4.1. Then,

χ 2 (t, y, k) = 1 γ ∂ 2 yy ϕ(t, y)X 2 (k).
Clearly III ǫ = O(1). Taking the limit in (4.9) we obtain (4.5), that ends the proof of the proposition.

Suppose that V (t, y) is a limiting point for some subsequence

V εn (t, y, k) = t 0 W ǫn (s, y, k)ds
where ε n → 0, in the sense described by Corollary 3.2. We can also assume that the respective sequence

( W ε (•)) ⋆-weakly converge to ρ(t, y) in L ∞ [0, +∞); L 2 w (R × T)
. For convenience sake we shall denote the subsequences by (V ε (•)), ( W ε (•)), respectively. We have V (t, y) = t 0 ρ(s, y)ds, for all t ≥ 0, and a.e. y. with D given by (4.6) and ρ 0 defined in (2.16).

Before showing the proof of the proposition we show how to finish, with its help, the proof of Theorem 2.2.

According to Proposition 4.3 V (t, y) is a weak solution of the heat equation ∂ t V (t, y) = D∂ 2 yy V (t, y) + ρ 0 (y) satisfying V (0, y) ≡ 0. According to part 3) of Corollary 3.2 we also have V (t, 0 ± ) ≡ 0, t ≥ 0. Hence,

V

2. 1 .

 1 Weak solution of the kinetic equation with an interface. In what follows we denote R -:= (-∞, 0). Consider W (t, y, k) := W (t, y, k) -T . It satisfies the equation (1.1) with the interface given by

  )

Proposition 4 . 2 .

 42 For any function ϕ ∈ C ∞ 0 (R + × R * ) we have+∞ 0 R ρ(t, y) ∂ t ϕ(t, y) + D∂ 2 yy ϕ(t, y) dtdy = 0. (4.5) with 0 < D := 1 γ T ω′ (k)X 1 (k)dk < +∞. (4.6)

Proposition 4 . 3 .

 43 For any function ϕ ∈ C ∞ 0 R+ × R * we have +∞ 0 RV (t, y) ϕ(t, y) + D∂ 2 yy ϕ(t, y) -ϕ(t, y)ρ 0 (y) dtdy = 0. (4.17

  )

  [START_REF] Basile | From a kinetic equation to a diffusion under an anomalous scaling[END_REF] we apply a version of the perturbed test function technique. Let ϕ ε (t, y, k) be determined by ϕ ε (t, y, k) = ϕ(t, y) + εχ 1 (t, y, k) + ε 2 χ 2 (t, y, k),

	y = 0.	(4.7)

  (t, y) = (y ′ )dy ′ , t, ±y > 0. (y ′ )dy ′ , t, ±y > 0,(4.18) which satisfies the conclusion of Theorem 2.2 for T = 0. The only thing yet to be shown is the proof of Proposition 4.3.Proof of Proposition 4.3. According to(3.8) we have∂ t V ε (t, y, k) + 1 ε ω′ (k)∂ y V ε (t, y, k) = γ ǫ 2 LV ǫ (t, y, k) + W ε (0, y, k) (4.19)and obviously V ǫ (0, y, k) = 0 a.e. Let ϕ ε (t, y, k) be given by (4.7). From (4.19) we can writeLV ε (t, y, k)ϕ ε (t, y, k) + W ε (0, y, k)ϕ ε y, k) dtdydk Lϕ ε (t, y, k) + W 0 (y, k)ϕ ε (t, y, k) dtdydk.Substituting from (4.7) we obtain that the term in the square brackets has the form (4.10), with I, II and III ǫ as given in (4.16). Taking the limit in (4.20) we obtain (4.7).

	+∞
	t ρ 0 This, combined with (4.16), implies that 0 1 √ 4πDs R± exp -(y -y ′ ) 2 4Ds -exp -(y + y ′ ) 2 4Ds ρ(t, y) = 1 √ 4πDt R± exp -(y -y ′ ) 2 4Dt -exp -(y + y ′ ) 2 4Dt 0 R×T ∂ t [V ε (t, y, k)ϕ ε (t, y, k)] dtdydk = +∞ 0 R×T V ε (t, y, k)∂ t ϕ ε (t, y, k) -1 ε ω′ (k)∂ y V ε (t, y, k)ϕ ε (t, y, k) + γ ε 2 = +∞ 0 R×T V ε (t, y, k) ∂ t ϕ ε (t, y, k) + 1 ε ω′ (k)∂ y ϕ ε (t, y, k) ρ 0 0 = γ + ε 2	(4.20)

T.K. acknowledges the support of the National Science Centre: NCN grant 2016/23/B/ST1/00492. SO's research is supported by ANR-15-CE40-0020-01 grant LSD. Both T.K. and S.O. were partially supported by the grant 346300 for IMPAN from the Simons Foundation and the matching 2015-2019 Polish MNiSW fund.


Comment on "Volume shrinkage of metal organic framework host induced by the dispersive attraction of guest gas molecules"

François-Xavier Coudert, a Alain H. Fuchs, a Alexander V. In a recent paper, 1 Joo et al. report from density functional theory calculations that H 2 adsorption in MOF-5 lead to a contraction of the metal-organic framework's unit cell. The authors contrast that finding with "the intuition that gas adsorption in a confined system (e.g., pores in a material) increases internal pressures and then leads to volumetric expansion", and describe the contraction as extraordinary. We would like to point out in this comment that adsorption-induced contraction of microporous materials is actually a very generic phenomenon, which has been demonstrated experimentally in a large range of microporous materials and explained theoretically. At the same time, we would like to acknowledge that the authors 1 should be credited for demonstrating this effect on the MOF-5 framework and characterizing its extent (of the order of 1% in volume). This is not atypical, if somewhat large, for inorganic microporous materials like zeolites, but small compared to metal-organic frameworks that show adsorption-induced structural transitions, such as the "breathing" MIL-53 (∆V/V ~ 40%).

Nonmonotonic deformation in the course of guest adsorption is the norm, and not the exception, in microporous materials: 2 these materials contract at low pressure (low loading), and expand at higher gas pressure (higher loading). This behaviour has been observed experimentally and documented dating back to 1940's 3 in zeolites, 4,5 microporous carbons, 6,7,[START_REF] Balzer | Proceedings of the Fifth Biot Conference on Poromechanics[END_REF] and metal-organic frameworks. [START_REF] Neimark | [END_REF] A good review of experimental measurements can be found in ref. 2. On the theoretical side, the contraction-expansion behaviour of microporous materials upon guest adsorption has also been well studied, both from the thermodynamic point of view based on the concept of adsorption stress 10,[START_REF] Neimark | [END_REF] and at the microscopic level using various molecular simulation tools. 11,12,13 This phenomenon is ubiquitous, and does not depend on the specifics of host-guest interactions, adsorbate distribution within the pores, nor does it require low temperature to observe. At low gas pressure, i.e. at low loading, contraction of the material is driven by the host-guest dispersion attractive interactions (provided the pores of the material are large enough to accommodate the guest). As the gas pressures increases, i.e. at higher loading, the adsorbed phase becomes denser and the host-guest distances diminish to the point that short-range repulsive interactions come to play and cause the material to expand in order to accommodate additional guest molecules. This is schematized on Fig. 1. contraction at low pressure, followed by expansion at higher pressure. This is here schematized using the structure of the highly compliant MIL-53 metal-organic framework.

  Neimark b a CNRS & Chimie ParisTech, 11 rue Pierre et Marie Curie, 75005 Paris, France b Department of Chemical and Biochemical Engineering, Rutgers University, Piscataway, New Jersey 08854, United States

Figure 1 .

 1 Figure 1. Nonmonotonic adsorption-induced deformation is typical in microporous materials:


Marta De Toni

Romain Jonchiere

Dr Pluton Pullumbi

François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Dr F.-X Coudert

Prof. Dr Alain H Fuchs

How Can a Hydrophobic MOF be Water-Unstable? Insight into the Hydration Mechanism of IRMOFs

Keywords: molecular dynamics, ab initio calculation, water, reaction mechanisms, metal-organic frameworks Introduction

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

is still relatively little data available on water adsorption in MOFs. 3 It is of great importance to understand the interaction of water with these materials, so as to better predict how they behave in real-life conditions, how they will age and keep (or loose) their surface area and adsorption properties over time. While some MOFs are stable in the presence of water, even when left in liquid water for long periods of time, others are unstable in the presence of water vapor or moist air. Such water-unstable MOFs include two archetypical materials, MOF-5 4 and HKUST-1, 5 both widely studied as candidates for gas separation and catalysis, but whose potential for use in industrial processes could be hindered by their limited hydrothermal stability. In the case of hydrophilic HKUST-1, it is easy to see that water binds strongly to unsaturated metal sites, 6 therefore leading to structural changes that ultimately lead to the framework's collapse; moreover, perfluorohexane deposition was recently demonstrated to improve the material's stability. 7 MOF-5, however, has no coordinatively unsaturated sites and presumably little attraction to water, making it harder to understand how it can be unstable in the presence of water vapor. In this paper, we study the interactions between water and an IRMOF framework and show that although the material itself is largely hydrophobic, it is not incompatible with a favorable hydration mechanism that leads to poor hydrothermal stability. We also show how ab initio (or first principles) molecular dynamics simulations can be used to investigate the stability of metal-organic frameworks upon exposure to water, and investigate their hydration mechanism.

The materials of the IRMOF family, discovered by Yaghi and co-workers, 4,8 feature oxide-centered Zn 4 O tetrahedral clusters as metal centers connected by dicarboxylate organic linkers and forming a 3D-connected cubic porous framework. IRMOF-1, also known as MOF-5, was the first known member of this family and is an archetypal metal-organic framework that has been widely studied and is now commercially available. Its organic linker is 1,4-benzenedicarboxylate (1,4-BDC), and it has a pore size of ~9 Å. Other members of this rich family feature different linkers, both smaller and larger in size, some with functional groups and heteroatoms. IRMOFs in general, and MOF-5 in particular, are known to be sensitive to moisture because it is believed that the weak zinc-oxygen coordination allows for an attack by water molecules. 9,10,11 Presence of water thus causes the framework to decompose, yielding lower surface area and poor adsorption properties. A number of recent publications have reported success in mitigating the issue by modifying the metal-organic framework with water-repellent functional groups, either by direct synthesis using functionalized linkers, 12,13 or by postsynthetic modification of the IRMOF framework. 14 Both routes allow one to tweak the hydrophobicity of the MOF, and yield materials that display almost no decline of their adsorption properties even after prolonged exposure to water.

In addition to these functionalized materials, some groups have tried to investigate the water-MOF interactions, both experimentally and theoretically. 3,6,15,16,17,18,19 In particular, an experimental study was reported by Sabo et al. on the conditions of hydrolysis of MOF-5 by water vapor, as well as its kinetics. 20 On the theoretical side, Greathouse et al. performed a constant-pressure, forcefield-based molecular dynamics of this same material to describe its interaction with water and the influence of water content on framework deformation and stability. 21 Finally, a combined experimental/theoretical study of the hydrothermal stability of various MOFs was published more recently by Low et al, 22 in which the authors conclude that framework stability in presence of water is mostly determined by metal-linker bond strength, which they measure using quantum mechanical calculations of MOFs cluster models. Their quantum chemistry calculations are based on the assumption that hydration takes place through linker displacement by a water molecule. In the present work, we investigate the interactions of water and the IRMOF framework in an unbiased and model-free approach, by using ab initio molecular dynamics. We previously highlighted the capabilities of this approach for the description of water adsorbed in hydrophobic and hydrophilic zeolites. 23,24 It combines the precision of quantum chemistry Density Functional Theory (DFT) calculations with a statistical sampling of phase space at non-zero temperature, allowing us to obtain more than minimum energy structures and transition states. We show here that it can thus be used to determine the hydration mechanism of a hybrid material. In addition to the ab initio MD, we have also performed forcefield-based Monte Carlo simulations, in order to ascertain the adsorption properties of IRMOFs.

Computational methods

Ab initio molecular dynamics of the anhydrous and hydrated MOF were performed with the Car-Parrinello approach 25 as implemented in the CPMD software package, 26 using the Kohn-Sham formulation of DFT. We employed the gradient-corrected BLYP functional. 27,28 A plane-wave basis set with an energy cutoff of 90 Ry was used and only the valence electrons were explicitly considered, with semilocal norm-conserving Troullier-Martins pseudopotentials taking into account the interactions between core and valence electrons. All molecular simulations were performed using periodic boundary conditions on a single cubic unit cell of the MOF. Hydrogen atoms were deuterated in order to improve decoupling of the dynamics of the ionic and electronic subsystems and allow for a larger time step in the integration of the equations of motion. A time step of 4 a.u. (~ 0.1 fs) and a fictitious electron mass of 400 a.u. were used.

For each water loading simulated, a starting configuration was extracted from the previously performed grand canonical Monte Carlo simulations (see below). A first equilibration period of at least 5 ps was performed on the system in the canonical ensemble, during which the temperature was controlled to 300 K by a Berendsen thermostat (with a time constant of 1000 u.a.). 29 Simulations were then run in the microcanonical NVE ensemble for respectively 185 ps (1 water molecule per unit cell) and 410 ps (4 H 2 O). The evolution of kinetic energy over time was checked to ensure that the NVE dynamics still corresponded to a temperature around 300 K.

In addition to ab initio MD simulations, we also performed forcefield-based Monte Carlo simulations in the grand canonical ensemble. The water molecules were described by the rigid, nonpolarizable TIP4P model. The IRMOF-0h structure was described using a forcefield derived from the one Snurr et al. developed for IRMOF-1, 30 where the atomic partial charges were slightly adapted to describe the acetylenedicarboxylate linker (instead of benzenedicarboxylate). These forcefield parameters are presented in the supplementary information. Lorentz-Berthelot rules were used to calculate MOF-water Lennard-Jones parameters. Series of GCMC simulations were performed at various values of water chemical potential, which was then related to water pressure by the same µ(P) relation as in ref. 24. All Monte Carlo runs were performed at 300 K on a 2x2x2 supercell to minimize the influence of periodic boundary conditions. Long-range electrostatic interactions were taken into account using the Ewald summation technique. In order to improve the efficiency of the calculations, electrostatic and repulsion-dispersion interaction energies between the rigid MOF framework and adsorbed water molecules were precomputed on a grid (with a grid mesh of 0.1 Å) and stored for use during the simulation.

Creation and characterization of the system

In order to characterize the hydrolysis mechanism of IRMOFs by ab initio molecular dynamics, we needed to find a member of the IRMOF family with a relatively small unit cell size so that current HPC resources would allow us to gather good statistics. IRMOF-1, also known as MOF-5, is the most studied member of the IRMOF family: its behavior in presence of water has been investigated, and functionalized MOF-5 structures have been devised that are relatively water-stable. However, with a cubic unit cell size of 25.8 Å containing 424 atoms, it is much too large for our ab initio molecular dynamics study. The only known IRMOF structure with a linker shorter than IRMOF-1 is Yaghi's IRMOF-0, 31 with formula (Zn 4 O)(acetylenedicarboxylate) 3 . However, IRMOF-0 experimentally forms as a doubly interpenetrating structure that exhibits no microporosity. We therefore created a hypothetical, non-interpenetrating analogue of IRMOF-0. We took as a starting point the experimental structure for one of the two interpenetrated networks, and ran a variable-cell DFT to minimize its energy (using the same functional, pseudopotentials and basis sets as described earlier for ab initio MD). The resulting structure is called IRMOF-0h. 32 It has a cubic spacegroup (𝑃43𝑚) and a unit cell parameter of a = 11.12 Å, close to the experimental value from the interpenetrated IRMOF-0 (a / 2 = 10.92 Å). It is presented in Fig. 1, and a CIF file of the structure is provided as part of the supplementary information. Grand Canonical Monte Carlo simulations of water adsorption in IRMOF-0h were then performed to characterize the affinity of water for this MOF. The adsorption-desorption isotherm cycle calculated by GCMC is shown in Fig. 2, with a change in horizontal scale at P sat = 5.2 kPa so that both the gas and liquid phase are clearly visible. This isotherm indicates that the ideal "classical" defect-free IRMOF-0h shows up as a hydrophobic material, with adsorption following a type V isotherm (zero slope around P = 0) and a saturation uptake in excess of 30 molecules per unit cell.

The hydrophobic nature of the solid is explained by the weak interactions between water and the carbon atoms of the framework, which are weaker than water-water interactions, as is the case of the archetypical water-graphene interaction. The step characteristic of the condensation transition happens at 20 MPa upon adsorption and that of drying at 2 kPa upon desorption, thus bracketing the water saturation pressure P sat with a large hysteresis loop characteristic of first-order condensationdrying transitions. Obviously, chemical events such as bond breaking cannot take place in these simulations based on classical forcefields. The isotherm shown in Fig. 2 thus describes what would happen in IRMOF-0h in absence of hydrolysis. We are then faced with an apparent contradiction. IRMOF-0h appears as a hydrophobic solid, while in real life the IRMOF materials are water unstable 8 . An examination of the water adsorption studies in several MOFs by Kaskel and coworkers 9 enables to understand that the hydrophilic/hydrophobic character of the inner surface of the MOF material should not be confused with the instability/resistance with respect to hydration. Indeed, these authors observe a Type I water adsorption isotherm in the case of HKUST-1, on account of the hydrophilicity of the free copper sites, accompanied by chemisorption and hydrolysis. In ZIF-8, a type V isotherm is observed, with a "remarkable chemical stability towards water". In the intermediate DUT-4 case, the adsorption branch of the isotherm is of type V, but the desorption branch shows that this aluminum based material, in which there are no free coordination sites available, experiences decomposition at high water loading. It turns out then that a hydrophobic MOF can be water unstable, a feature that was overlooked in some previous studies. 6,7,8 The following ab initio simulations provide some insight into this phenomenon.

Ab initio molecular dynamics at low water loading

A first ab initio simulation was run with a single water molecule adsorbed in a unit cell of IRMOF-0h, for a total of 185 ps of NVE molecular dynamics. For most of the trajectory's length, the water molecule diffuses rapidly in the nanopores of the material. However, there are also periods of time during which it gets closer to one zinc ion (down to 2 Å, see fig. 3). The histogram of Zn-O distances (fig. 4) shows a clear peak at about 2.2 Å. This state corresponds to the water molecule being bound to Zn 2+ ion, i.e. entering its coordination complex. To further characterize this state during the simulation, we decided on a simple geometric threshold at the minimum in the distribution of distances: we consider that d(O-Zn) ≤ 2.8 Å indicates a bound state, while larger distances correspond to a free water molecule. With this criterion, we see that over the 160 ps of the dynamics, three events of water attack on the metal centre occurred. Each of the bound states produced (see Table 1) are relatively short-lived, lasting no more than 3.0 ps (see Table 2). Hence, water attack on the zinc cation is a relatively frequent event, is reversible, and the free ⇄ bound equilibrium strongly favors the free state at this low water loading. The existence of this bound state is consistent with the hydrophobic nature of the host matrix and its type V isotherm.

Indeed, hydrophobicity does not mean there are no interactions (or purely repulsive interactions), but only that host-guest interactions are weaker than guest-guest interactions. 33 Hence it is possible for a bound state to exist, but it must be thermodynamically unfavourable compared to the freely diffusing molecule for the host to be considered hydrophobic.

In order to quantify this bound state further, we computed the potential of mean force as a function of Zn-O water distance from our dynamics. While this certainly bears some uncertainty due to the relatively small statistical sampling compared to the frequency of the reaction, it does give an estimate of the free energy difference between the two states, and the barrier height. The calculated potential of mean force is plotted on Fig. 4. It shows a free energy barrier for the attack of ∆F ≠ = 3.2 kT (~8 kJ/mol) and a free energy difference between metastable states of 1.3 kT (~3.2 kJ/mol) in favor of the water-free state.

Furthermore, we characterized this state of one water molecule bound to the Zn 2+ ion. The cation has a tetrahedral environment in the nonhydrated material (state A in the overall sketch on Fig. 5), with 4 oxygen atoms coordinated to it: 3 from the carboxylate linkers, and 1 being the central O atom of the Zn 4 O cluster. When the water molecule binds to the Zn 2+ , none of these atoms exit the coordination shell, and the Zn 2+ thus has a trigonal bipyramidal environment (state B in Fig. 5).

This is rather typical of Zn(II), which is commonly encountered in high-coordination structures, including pentacoordinate states. 34 , 35 This structure is moreover similar in geometry to the "hydrated" state hypothesized by Low et al. 22 as the initial state of their ligand displacement mechanism. We further confirmed that this state is indeed a stable state, and obtained its geometric characteristics (given in Table 1) by an energy minimization with a variable cell (using the same DFT parameters as for the MD run). 1. Geometrical characteristics of Zn 2+ environment each state described in Fig. 5.

Effect of higher water loading

In order to investigate the effect of hydration level on the stability of IRMOFs, we ran ab initio MD simulations with a loading of 4 water molecules per unit cell of IRMOF-0h. The attacks of water molecules on the Zn 4 O metal cluster observed for 1 and 4 water molecules, and the lifetimes of the bound states produced, are presented in Table 2. It is clear that the presence of additional water molecules in the IRMOF-0h unit cell stabilizes the bound state and displaces the A ⇄ B equilibrium towards water-bound state B, precursor of hydrolysis. For the 4 water molecules system, which ran for a total of 410 ps, a single water attack event was observed at 34.5 ps, after which the bound state B was never unformed. Because of this single event, we could not estimate the free energy difference between the A and B states at a water loading of 4 H 2 O / u.c., but only conclude that B is much more stable than A in these conditions.

To analyse the reasons for this stabilization, we characterized the hydrogen bonds network of the adsorbed water. We used a simple geometric criterion: 36 Besides the bipyramidal water-bound state observed upon H 2 O attack of the metal center, the ab initio MD for 4 adsorbed water molecules revealed a second hydrated state of the IRMOF-0h. In this state, labeled C and represented in Fig. 6, one of the carboxylate linker is no longer bound to the Zn 2+ cation, which thus regains its tetrahedral coordination (water molecule, Zn 4 O central oxygen and 2 carboxylate oxygens). Because this state has a carboxylate linker no longer coordinated to two metal cations, i.e. with a "dangling oxygen", it clearly constitutes a second step in the hydration mechanism (Fig. 5), i.e. the displacement of the linker. The C state is short-lived and the second reaction (B ⇄ C) is again reversible (see Table 2), and in fact in rapid equilibrium between what can be seen as two labile states. Moreover, we observe from the evolution of the various Zn-O carboxylate distances pictured over time on Fig. S3 that the linker displaced by the water molecule is not the same in every attack. This hydrated state can thus be described as a dynamic

equilibrium B ⇄ C 1 ⇄ B ⇄ C 2 ⇄ B ⇄ C 3 .
We expect that this situation would be somewhat modified in a real experimental system of macroscopic size: in constant-pressure conditions, the deformation of the unit cell gives the linker full latitude to move, and if a system larger than a single unit cell were considered, more than one metal cluster could be attacked by water. As a consequence, the linker once displaced would later move away from the zinc structure, eventually leading to the collapse of the structure. (the water and central oxygen atoms being on apical positions), while in the almost tetrahedral state C this angle is closer to 109.4°. The system can clearly be seen to switch between these two welldefined states (Fig. S3). The potential of mean force along this reaction coordinate is featured in Figure 7, which confirms the small free energy barrier between the two states. It also indicates that under these conditions, structure C with the displaced ligand is slightly favored over hydrated state B.

Overall mechanism and comparison with previous studies

The overall mechanism that we derived from our ab initio MD simulations is presented above in Fig. 5. It is important to note that this mechanism is not an hypothesis of our work, but was derived entirely from direct first-principles simulations at various water loading. Hydration proceeds in two successive steps: (i) water binding to the metallic cluster; (ii) linker displacement, where one water molecule actually replaces one of the linker arms. These two steps take place very rapidly, since they can be observed on the timescale of ab initio MD, provided that there is enough water present in the porous framework (here 4 molecules per unit cell corresponds to 12 weight%). One water molecule binds to the metallic cluster and the extra water molecules stabilize both the hydrated and the linker-displaced states so that the energy barrier for the hydration reactions become of the order or less than kT. Thus, we showed that the hydration and linker displacement mechanisms cannot be well-described by cluster calculations with a single water molecules: a full unit cell has to be used, O linker . This contrasts with our results presented above, where we characterized two possible hydrated states at low to medium water loading, none of which displays a rupture (or even an marked elongation) of the Zn-O central bond. In fact, this bond appears to be the stronger of the two coordinative bonds, and the Zn 2+ cation will prefer keeping a coordination number of five, with a trigonal bipyramid geometry, rather than break the Zn-O central bond. We believe this difference with Greathouse's findings is due to a bias of the classical forcefield, which has a tendency to describe the Zn 2+ as tetrahedrally coordinated.

Conclusions

The inner surface of porous IRMOFs, modelled in this work by the smaller IRMOF-0h, mainly exhibit non polar functional groups together with Zn ions that are coordinated to oxygen atoms.

This surface has an overall hydrophobic character, and this explains the type V water adsorption isotherm observed in a classical Monte Carlo simulation. The Zn 4 O clusters act as "weak hydrophilic defects", as defined in ref. 37 and 38. When only one water molecule is present in the pore, the bound state is frequently visited but it is less stable than the water free state, and no chemisorption of any kind is taking place. At a higher water loading, a water cluster can be formed at the Zn cluster site and this was shown to stabilize the water bound state, which very rapidly transforms into a linker-displaced state, where water has fully displaced one arm of a linker and which corresponds to the loss of the material's fully-ordered structure. Thus an overall hydrophobic MOF material can also become water unstable, a feature that was not fully understood until now.

In future works, we aim at extending our ab initio molecular dynamics simulations to the isothermal-isobaric (N,P,T) ensemble with a varying unit cell, so that the complete hydration event could be simulated. This kind of study will be extended to other families of metal-organic frameworks, in order to evaluate the influence of metal-linker coordination chemistry on the stability of these materials. We will also investigate the stabilization of functionalized IRMOF variants such as CH 3 -MOF-5, for which the steric effect of the methyl group is not enough to justify the increased hydrothermal stability of the material.

Å

Radial Distribution

Figure 1 .

 1 Figure 1. Representation of the cubic structure of IRMOF-0h. Carbon atoms are in black, oxygen in red, Zn-centered tetrahedra in blue.

33

 33

Figure 2 .

 2 Figure 2. Adsorption-desorption isotherms of water in IRMOF-0h calculated from GrandCanonical Monte Carlo simulations at 300 K.

Figure 3 .

 3 Figure 3. Distance between the oxygen of the water molecule and the closest Zn 2+ ion in a simulation of one H 2 O in IRMOF-0h at 300 K. The threshold of 2.8 Å chosen to characterize the bound state of the water molecule is indicated as a red dashed line. Each occurrence of a bound state during the dynamics is indicated by a thick horizontal line, labelled with the index of the Zn 2+ ion involved.

Figure 4 .

 4 Figure 4. Histogram of water-Zn distance for one water molecule adsorbed in IRMOF-0h at 300 K (upper panel) and potential of mean force (lower panel).

Figure 5 .

 5 Figure 5. Different states of a Zn 2+ coordination sphere of IRMOF-0h with zero or one bound water molecules. See text for details on each state.

Figure 6 .

 6 Figure 6. Snapshot of a water molecule bound to the Zn 2+ ion of IRMOF-0h, with tetrahedral coordination (hydrated state C). The other 3 water molecules which form a hydrogen-bonded network with the first one (see text) are not shown.

Figure 7 .

 7 Figure 7. Potential of mean force as a function of O c ZnO water angle in the MD simulation of 4 H 2 O molecules in IRMOF-0h at 300 K, showing the equilibrium between states B and C.

  and the influence of the material's water loading need to be taken into account. The studies of ref 22 for instance were performed with a single water molecule and gave a very high value of the energy barrier for linker displacement (∆E ≠ ~ 48.6 kJ/mol). Finally, the comparison of the present work with the results of previous theoretical work by Greathouse et al 21 , who performed constant pressure forcefield-based classical MD of hydrated IRMOF-1, reveals stark differences in the structures obtained for the bound Zn 2+ -water states. From the results of their forcefield-based simulations, Greathouse et al concluded that water displaced rather indifferently two types of oxygen atoms initially coordinated to the Zn 2+ ions: the central atom of the Zn 4 O tetrahedron, which we call O central , and those of the linkers' carboxylate groups,

  Figure 2

  Figure 5

Table

  coordination sphere of IRMOF-0h with zero or one bound water molecules. See text for details on each state.

		A	B	C
	geometry	tetrahedral	trigonal bipyramid	tetrahedral
	d(Zn-O central )	2.00 Å	1.98 Å	2.04 Å
	d(Zn-O linker )	2.00 Å	2.04 Å	1.97 Å
	d(Zn-O water )	-	2.06 Å	2.17 Å

4 water molecules (total: 410 ps)

  d OO ≤ 3.5 Å, θ (OOH) ≤ 40°. During all the time spent with a water molecule bound to Zn 2+ (state B), the other three water molecules in the system tend to stay in the vicinity of the bound molecule, forming a weak hydrogen bonds network.Over the course of the entire simulation (276 ps), this hydrogen bonds network reorganizes several times, leading to various geometries featured in Fig.S2. However, we observed at least one water molecule hydrogen-bonded to the bound one at all times. The existence of this network of water molecules accounts for the stabilization of the bound Zn-H 2 O state with increasing water loading.

	1 water molecule (total: 185 ps)
	H 2 O index	Zn index	Attack time	Departure time	Lifetime of bound state
	H 2 O(1)	Zn(2)	47.5 ps	50.5 ps	3.0 ps
	H 2 O(1)	Zn(4)	105.0 ps	105.3 ps	0.3 ps
	H 2 O(1)	Zn(1)	117.0 ps	118.4 ps	1.4 ps
	H 2 O index	Zn index	Attack time	Departure time	Lifetime of bound state
			34.5 ps	_	≥ 375 ps
			State B (trigonal
			bipyramid)
			34.5 ps	97 ps	62.5 ps
			206 ps	222 ps	16.0 ps
			264 ps	282 ps	18.0 ps
			311 ps	338 ps	27.0 ps
			State C
			(tetrahedral)
			97 ps	206 ps	109 ps
			222 ps	264 ps	42 ps
			282 ps	311 ps	29 ps
			338 ps	_	≥ 72 ps

Table 2 .

 2 Analysis of the attacks of water molecules on Zn 2+ ions, and the lifetime of the zinc-water bound states.

Acknowledgments

The authors thank Rodolphe Vuilleumier and Anne Boutin for their technical expertise and many fruitful discussions. This work was supported by the Agence Nationale de la Recherche under contract ANR-2010-BLAN-0822-02, and performed using HPC resources from GENCI-IDRIS (project 096114).