Sami Boutamine
email: sami.boutamine@utc.fr

Dan Istrate
email: dan.istrate@utc.fr

Jérôme Boudy
email: jerome.boudy@telecom-sudparis.eu

Smart Sound Sensor for distress situation and number of presented person

Keywords: Speaker diarization, home monitoring, sound recognition, aal, aha, signal processing

In order to allow older adults to active and healthy ageing with comfort and security we have already developed a smart audio sensor being able to recognize everyday life sounds in order to detect activities of daily living (ADL) and distress situations. In this paper, we propose to add a new functionality by analyzing the speech flow in order to detect the number of person in a room. The proposed algorithms are based on speaker diarization methods. This information is useful in order to better detect activities of daily life but also to know when the person is home alone. This functionality can also offer more comfort through light, heating and air conditioning adaptation to the number of persons.

I. INTRODUCTION

Nowadays, technologies at the service of intelligent spaces are constantly developing and interacting with everyday life objects. They exploit video, audio signals and environmental data to locate people, to recognize their gestures, to interact with the people in order to offer comfort and to help people with disabilities. An important domain of study is the Ambient Assisting Living (AAL) which try using the new technologies to compensate the disabilities related to the age (visual, hearing, cognitive).

We have already developed a smart audio sensor [START_REF] Robin | Remote monitoring, distress detection by slightest invasive systems: Sound recognition based on hierarchical i-vectors[END_REF] allowing the detection of distress situations but also of people activity through sound environment analysis. This proposed sensor analyse the sound environment and is able using different techniques to recognize 18 sound classes in a continuous audio flow.

In this paper, we propose a new functionality to the already presented sensor by adding the possibility to detect the number of persons present in a room through speech analysis. In fact, knowing the number of persons help to know if the elderly people is home alone and to adapt the distress situations identification. Also, knowing the presence of other persons allow to follow the intervention at home of different services (catering, cleaning services, …). Otherwise, combining the speech with sound analysis allow a better activity daily life identification in order to allow an Active and Healthy Ageing (AHA) of young older adults.

Additionally, this new functionality allow also offering comfort and energy consumption reduction by adapting the light quantity, the heating or the air conditioning systems. This work is a part of the multi-partner national project CoCAPs (FUI type). In this article, we present the work and results obtained from the detection of the number of speakers where several tools and methods of sound processing are used including speaker diarization. This paper is organized as follows, we first present the speaker diarization system, then, we highlight the developed application (detection of the number of speakers) by detailing the implementation of the application of sound detection system (speech/no-speech) through the use of functionalities offered by the LIUM_SpkDiarization toolkit. Finally, we present and discuss the first test results of the developed application.

II. SPEAKER DIARIZATION SYSTEM

Speaker diarization is to determine "who spoke when?" in an audio record that contain speech, music or noise segments. The signal audio is splitted into homogeneous speech segments, according to the speaker identity with no prior knowledge about it.

The principal modules of diarization system are composed of parametrization, speech activity detection or voice activity detection, speaker segmentation, speaker clustering and resegmentation.

Already available tools are LIUM_SpkDiarization [START_REF] Meigner | An Open Source Toolkit For Diarization Sylvain Meignier[END_REF], audioSeg, DiarTK, and SHoUT for this purpose.

A. Parametrization

Mel-Frequency Cepstral Coefficients (MFCC), Perceptual Linear Prediction coefficients (PLP) and Linear Frequency Cepstral Coefficients (LFCC) are the most common features used to extract the most important information from the signal, sometimes used with their first and/or second derivatives. Ideal features should have some important properties:

• They have to emphasize the difference between classes (class separability). • Their intra-class variances must be minimal.

• They have to be robust to noise trouble, conserving the class separability as far as possible. • The total number of detected features (i.e. features density) should be sufficiently large to reflect the frames content in a compact form. • A high correlation between features should be avoided, as much as possible.

The MFCCs are widely used in automatic speech and speaker recognition, and used for this current work.

B. Speech activity detection

The speech activity detection plays a very important role in the whole diarization process, it is performed to separate speech, no-speech and silent frames. The output of other sub tasks highly depends on the precision of this task. The goal is to keep only relevant information for speech modeling and to reduce the amount of computation needed for the following treatments. Many approaches allow the detection of speech, segmentation using Hidden Markov Models (HMMs) [START_REF] Renevey | Entropy based voice activity detection in very noisy conditions[END_REF], Deep Neural Networks (DNN) [START_REF] Ryant | Speech activity detection on youtube using deep neural networks[END_REF]. In this work, we used an algorithm based on the Wavelet Transform [START_REF] Istrate | Medical Telemonitoring System Based on Sound Detection and Classification[END_REF].

C. Speaker segmentation

The aim of the speaker segmentation is to find speaker change points in a given speech signal using the symmetric Kullback Leibler distance (KL2), the generalized likelihood ratio (GLR) or the Bayesian information criterion (BIC) distance computed using Gaussians with full covariance matrices. In such systems, the speech signals are windowed for a short duration of 25 -30ms.

In such systems like LIUM_SpkDiarization toolkit [6], the segmentation is done in two steps. A first pass on the signal is performed to detect breaks (changes in speakers), using the GLR (Generalized Likelihood Ratio) measure defined in "Eq. ( 1)". The GLR measure introduced by Gish [START_REF] Gish | Segregation of speakers for speech recognition and speaker identification[END_REF] is a likelihood ratio between two hypotheses H 0 and H 1 .

H 0 : The two sequences x i and x j are produced by the same speaker x, which, in this case, the model M(µ, ) corresponding to x = x i ∪ x j would allow a better representation of x i and x j . H 1 : The two sequences x i and x j are produced by two different speakers, which case the two model M i (µ , ) and M j (µ , ) would be better suited to represent x i and x j .

The likelihood test is thus formulated by the ratio of the two hypotheses:

, = ( , (µ, Σ)) ( , (µ , Σ )) ( , M µ , Σ ) (1)
Where ( , (µ, ))corresponds to the likelihood of the sequence x = x i ∪ x j given the model M(µ, ), and( , (µ , )) ( , (µ , )) the likelihood that the x i and x j were produced by two different speakers.

A second pass, allows to refine the segmentation obtained during the first pass by grouping consecutive segments that maximize a likelihood score using a discriminating measure BIC (Bayesian Information Criterion) defined in "Eq. ( 2)" below. BIC is a metric highly appreciated for its simplicity and efficiency.

= -2 ln( ) + ln( )

With L the likelihood of the estimated model, N the number of observations in the sample and k the number of free parameters of the model.

D. Speaker clustering

Clustering is done using an unsupervised method called hierarchical agglomerative clustering (HAC). The goal of speaker clustering is to associate segments from an identical speaker together. Speaker clustering ideally produces one cluster for each speaker with all segments from a given speaker in a single cluster.

The initial set of clusters is composed of one segment per cluster. Each cluster is modeled by a Gaussian with a full covariance matrix. ∆BIC measure is employed to select the candidate clusters to group as well as to stop the merging process. The two closest clusters i and j are merged at each iteration until ∆BIC i,j > 0. ∆BIC is defined in "Eq. ( 3)" [6].

∆ , = " # $" % & log|Σ| - " # & log|Σ | - " % & log*Σ * -λ + (3) + = , & -+ . (.$,) & + log / + / (4
) Where |Σ i |, |Σ j | and |Σ| are the determinants of gaussians associated to the clusters i, j and i + j. λ is a parameter to set up. The penalty factor P "Eq. ( 4)" depends on d, the dimension of the features, as well as on n i and n j , referring to the total length of cluster i and cluster j respectively.

This penalty factor only takes the length of the two candidate clusters into account where as the standard factor uses the length of the whole data.

E. Re-segmentation

Re-segmentation is the final stage of the process, in which the rough boundaries of diarization systems that rely on segment clustering of an initial uniform segmentation are refined based on a frame-level. The most common approach is a Viterbi re-segmentation with MFCC features [START_REF] Meigner | An Open Source Toolkit For Diarization Sylvain Meignier[END_REF].

III. APPLICATION

The purpose of the application developed is the detection of the number of speakers in a room, a meeting room or an office, the development was realized on the Raspberry Pi3 Model B board.

In order to achieve our objective, we used the speaker diarization technique, whose goal is to segment the audio signal into small homogeneous regions containing only speech and which belongs to one and only one speaker. The number of speakers corresponds to the number of segments group obtained from each speaker.

A. Application architecture

The basic architecture of the application is shown in (Fig. 1).

Figure 1. Application Architecture

It is an application that works in parallel, using two threads. In the first thread an application proceeds in a continuous way to capture the sound signal and differentiate it between the speech and the other sounds of everyday life whose objective is the recording of the segments containing only speech.

In the second thread, the program gets back and groups the segments of saved speech for a predefined duration in order to move to the segmentation phase using the LIUM_SpkDiarization toolkit. The output of this phase is represented by a file containing groups of segments with the speech of each speaker. The number of groups corresponds to the number of people detected automatically by speaker diarization.

B. Sound detection system (speech/no-speech)

After analyzing the first results of the toolkit LIUM_SpkDiarization applied to a complete audio signal containing speech and human sounds of everyday life, it has been noticed that sounds, which are different from the speech, have a negative influence on the results. Indeed a sound detection application was applied before using the toolkit LIUM_SpkDiarization, whose purpose is to filter the audio signal by removing all other sounds different from the speech.

The sound detection aims to detect, and separate speech events from other sound events in the continuous audio flow.

The classification of sound/speech is based on two Gaussians Mixtures Models (GMM). One class for speech and another one for the sounds of everyday life.

C. LIUM_SpkDiarization toolkit

LIUM_SpkDiarization is an open source diarization toolkit designed for extraction of speaker identity from audio records with no information before about the analysed data (number of speakers, etc.). LIUM could identify speaker's speech segments at excellent level. LIUM_SpkDiarization was developed by LIUM (Computer Science Laboratory of Le Mans) for the French ESTER2 evaluation campaign, the architecture of this system is shown in (Fig. 2). LIUM_SpkDiarization is written in Java, to minimize dependency problems with the various operating systems and libraries. This tool is distributed as one selfcontained JAR archive, which can be run directly [6].

LIUM_SpkDiarization comprises a full set of tools to create a complete system for speaker diarization, going from the audio signal to speaker clustering based on the CLR/NCLR metrics. These tools include MFCC computation, speech/nospeech detection, and speaker diarization methods [START_REF] Eva | Comparison of Diarization Tools for Building Speaker Database[END_REF].

1) System description

The diarization system provided by LIUM starts its processing by the computation of the acoustic parameters (MFCC) which are compute directly from the signal by the program using Sphinx4 (one of the included third-party packages) [6]. The features are composed of 13 MFCCs with coefficient C0 as energy for segmentation based on BIC (Bayesian Information Criterion), BIC clustering and segmentation based on Viterbi decoding. The features are reduced to 12 MFCCs completed by ∆ coefficients (coefficient C0 is removed) for speech detection, Gender and band width detection.

The speaker's segmentation is composed of two phases based on GLR (Generalized Likelihood Ratio) that identifies the instantaneous change points and BIC (Bayesian Information Criterion) distance metric for the speaker consecutive segments fusion. The speakers are modelled by Gaussian distribution with full covariance matrix in the segmentation and clustering phases.

LIUM_SpkDiarization system finally performs another HAC (hierarchical agglomerative clustering) using normalized Cross-Likelihood Ratio (CLR) or Integer Linear Programming (ILP) proposed, where i-vectors were used to model and measure the similarity between clusters [START_REF] Rouvier | An Open-source State-of-the-art Toolbox for Broadcast News Diarization[END_REF].

IV. EXPERIMENTS AND RESULTS

The first tests done to evaluate the developed application "detection of the number of persons" are performed in offices, taking into account the number of speakers and the duration of the audio segments to be processed.

The properties of the processed audio signal are given in (Table I). The Table II, presents the performance of the sound detection system (speech / no-speech), using an audio signal lasting 7min containing speech and sounds. The results obtained are shown by the curves in (Fig. 3) and in (Table III), the numbers represent the percentages of the algorithm performance in relation to the number of speakers and the durations of the audio segments.

The results represent the performance of the application "detection of the number of persons", applying the sound classification speech / no-speech to remove any different sounds of speech, followed by the speaker diarization tool (LIUM_SpkDiarization), the suppression of sound has improved the results. From the first results, we notice that the algorithm performance is better for the longer durations (10 minutes).

The variation of voice attenuation for the same speaker present in the studied signal has a negative influence on the obtained results.

The tables (Table IV), (Table V) and (Table VI) represent the tests results performed for a signal audio lasting 3 minutes, 5 minutes and 10 minutes respectively, containing three speakers.

S1(F) S2(F) S3(F) Speaker1(F) 1min-03sec X X Speaker2(F) X 59sec X Speaker3(F) X X 48sec TABLE V. EXAMPLE TESTS (5 MINUTES, 1 WOMAN, 2 MEN) S1(M) S2(M) S3(F) Speaker1(M) 51sec X X Speaker2(M) X 1min-58sec X Speaker3(F) X X 1min-29sec
S1(M) S2(F) S3(F) Speaker1(M) 4min-23sec X X Speaker2(F) X 2min-02sec X Speaker3(F) X X 2min-16sec
S1, S2 and S3 correspond to the speaker label found automatically by the algorithm, M and F mean respectively Male and Female, the duration represent the time that each speaker has spoken.

The algorithm can make the difference between man and woman with an error rate equal to 0%, and this represents a very important information in our application, especially for the recognition of person in a room.

The evaluation of the application is underway, notably for a larger number of speakers and in different environments such as living rooms and meeting rooms, the aim of which is to improve the results obtained.

V. CONCLUSION AND PERSPECTIVES

This paper proposes a new functionality to an existing smart audio sensor being able to recognize everyday life sound for ADL and distress detection. The new functionality allow the estimation of the number of speaker in the speech flow, information useful for ADL detection but also for adaptation of distress detection system when the person is home alone.

The results of this work corresponds to the tests of the speaker number detection application based on both speaker diarization (LIUM_SpkDiarization tool) and the application of the sound classification (Speech / Sound), by removing any sounds different from speech before applying the speaker diarization tool, whose goal is to improve the results obtained by LIUM_SpkDiarization.

For the future, the results will be improved by working on the sensor location optimization and the audio signal filtering. The study of cases where several speakers speak at the same time is envisaged. The use of speaker recognition methods is potentially considered as to improve the remote monitoring for elderly people in a room.

Figure 2 .

 2 Figure 2. Architecture LIUM_SpkDiarization [6]

TABLE I .

 I SPEECH CORPUS

	Language	French
	Sampling rate	16 Khz
	N° of channels	1 (16-bits mono channel)
	Speech domain	Conversational speech

TABLE II .

 II SOUND DETECTION SYSTEM PERFORMANCE

		Speech	Sound
	Speech	89,66%	10,34%
	Sound	16,67%	83,33%

TABLE III

 III

		.	TEST RESULTS
		3min	5min	10min
	1 Speaker	70%	90%	90%
	2 Speakers	50%	70%	90%
	3 Speakers	90%	70%	80%
	Total	70%	76.66%	86.66%
	Diarization performance	77.77%
	Figure 3. Application performance

TABLE IV

 IV

	.	EXAMPLE TESTS (3 MINUTES, 3 WOMEN)

ACKNOWLEDGMENTS

The authors would like to thank BPI France, the Regional Councils of Limousin and Rhône-Alpes associated with the ERDF program, the departmental council of Isère, and the Bourges agglomeration community for their financial support to the project CoCAPs.

The CoCAPs project, from FUI N ° 20, is also supported by the poles of competitiveness S2E2 and Minalogic.


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

Romain Gaillac

Pluton Pullumbi

Melting of Zeolitic

  L'archive ouverte pluridisciplinaire

Introduction

Metal-organic frameworks (MOFs) appeared almost 30 years ago [START_REF] Hoskins | Infinite polymeric frameworks consisting of three dimensionally linked rod-like segments[END_REF] in the literature, and much of the research effort and publications have been focused on the study of crystalline structures [START_REF] Furukawa | The chemistry and applications of metal-organic frameworks[END_REF] . These microporous or mesoporous crystalline architectures, in addition to being beautiful chemical objects, exhibit very large specific surface areas and pore volumes, making them good candidates for applications such as gas sorption or heterogeneous catalysis. Moreover, their crystallinity makes them relatively easy to characterize and allows the determination of their structure, typically by X-ray crystallography. Nonetheless, there has been a recent shift, with increasingly more examples of MOFs that show defects and disorder, and the realization that this is not always detrimental to their physical or chemical properties -and thus, to their performance in potential applications.

There are even studies that focus on "defective by design" strategies to explore further the configurational space of organometallic structures [START_REF] Bennett | Interplay between defects, disorder and flexibility in metal-organic frameworks[END_REF] . In particular, several amorphous MOFs have been studied in the last few years and proposed as candidates for applications such as controlled release, irreversible trapping of harmful substances, or creating new optically active glass-like materials [START_REF] Bennett | Amorphous metal-organic frameworks[END_REF] .

There are several ways to obtain amorphous MOFs, ranging from direct synthesis [START_REF] Zhao | Nanoporous transparent MOF glasses with accessible internal surface[END_REF] to the amorphization of crystalline MOFs by melting and quenching, [START_REF] Bennett | Hybrid glasses from strong and fragile metal-organic framework liquids[END_REF] application of mechanical pressure, [START_REF] Coudert | Responsive metal-organic frameworks and framework materials: under pressure, taking the heat, in the spotlight, with friends[END_REF] or ball-milling [START_REF] Bennett | Ball-milling-induced amorphization of zeolitic imidazolate frameworks (ZIFs) for the irreversible trapping of iodine[END_REF] .

It has been known for 15 years that inorganic zeolites can undergo temperature-induced amorphization. [START_REF] Greaves | The rheology of collapsing zeolites amorphized by temperature and pressure[END_REF][START_REF] Colligan | Synchrotron X-ray powder diffraction and computational investigation of purely siliceous zeolite Y under pressure[END_REF][START_REF] Haines | Topologically ordered amorphous silica obtained from the collapsed siliceous zeolite, silicalite-1-F: a step toward "perfect" glasses[END_REF] Zeolitic Imidazolate Frameworks (ZIFs) are a subgroup of particularly thermally and chemically stable MOFs exhibiting zeolitic topologies. They are formed by assembly of tetrahedra based on a metal ion and four imidazolate ligands (possibly substituted). In 2010 was published the first report of an amorphous zeolitic imidazolate framework formed by heating Zn-based framework ZIF-4. [START_REF] Bennett | Structure and properties of an amorphous metal-organic framework[END_REF] Then Zn-based ZIF-1, ZIF-3 and Co-based ZIF-4 were shown to amorphize upon temperature increase, while ZIF-8, ZIF-9, ZIF-11, ZIF-14 and ZIF-β qtz would decompose without undergoing a phase transition. [START_REF] Bennett | Thermal amorphization of zeolitic imidazolate frameworks[END_REF] The only common point to the latter are their substituted ligands compared to the former which are composed of non-substituted imidazolate cycles. In Ref. 6 Bennett et al. detail the differences in behaviour between ZIF-8 and ZIF-4 in terms of calorimetric and of the ligand and the topology on the melting temperature for ZIFs which can undergo thermal amorphization. [START_REF] Bennett | Melt-quenched glasses of metal-organic frameworks[END_REF] For instance, ZIF-4 was compared to TIF-4 and ZIF-62 which possess the same topology but have slightly larger ligands. Both, TIF-4 and ZIF-62, then have shown lower melting temperatures than ZIF-4.

Recently, we studied the melting mechanism of ZIF-4, [START_REF] Gaillac | Liquid metal-organic frameworks[END_REF] using a combined experimental and computational study, and in particular relying on in situ variable temperature X-ray, ex situ neutron pair distribution function experiments, and first-principles molecular dynamics simulations. Here, we present a theoretical study on the melting mechanism of three different ZIF frameworks, namely ZIF-8, ZIF-4 and ZIF-zni, focusing on the influence of the topology of the crystalline MOF on the melting behaviour and the structure of the liquid phase. From the three frameworks chosen, two are known to undergo thermal amorphization (ZIF-4, ZIF-zni), while the other one collapses before melting (ZIF-8). ZIF-8 is composed of Zn(mIm) 2 tetrahedra linked by Zn-N coordination bonds and has the same topology as the zeolite sodalite (SOD) [START_REF] Coudert | Molecular mechanism of swing effect in zeolitic imidazolate framework ZIF-8: continuous deformation upon adsorption[END_REF] . ZIF-4 and ZIF-zni are chemically very similar to ZIF-8, except for Zn(Im) 2 tetrahedra instead of Zn(mIm) 2 . ZIF-4 is porous and shares its topology with the mineral variscite (cag topology) while ZIF-zni has a non-porous dense zni topology [START_REF] Tan | Chemical structure, network topology, and porosity effects on the mechanical properties of zeolitic imidazolate frameworks[END_REF] . ZIF-8, because of its SOD topology, is a very open and porous framework.

Methods

First principles molecular dynamics

The behaviour of zeolitic imidazolate frameworks as a function of temperature was studied by means of density functional theory (DFT)-based molecular dynamics (MD) simulations, using the Quickstep module [START_REF] Vandevondele | Quickstep: fast and accurate density functional calculations using a mixed Gaussian and plane waves approach[END_REF] of the CP2K software package. 19 We used the hybrid Gaussian and plane wave method GPW [START_REF] Vandevondele | Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach[END_REF] as implemented in CP2K. The simulations were performed in the constantvolume (N,V, T ) ensemble with fixed size and shape of the unit cell. A timestep of 0.5 fs was used in the MD runs, the temperature was controlled by velocity rescaling. [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] The unsually large temperature at which the simulations were performed requires careful fine-tuning of the simulation protocol. In particular, the exchange-correlation energy was evaluated in the PBE approximation, [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and the dispersion interactions were treated at the DFT-D3 level. [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 Quickstep module uses a multi-grid system to map the basis functions onto. We kept the default number of 4 different grids but chose a relatively high plane-wave cutoff for the electronic density to be 600 Ry, as already used in Ref. 24

Trajectory analysis

The Lindemann ratio ∆ is computed from the width of the first peaks in the different partial radial distribution functions as a measure of the fluctuation of atomic positions and interparticular distances:

∆ = FWHM d 0 ( 1
)
where FWHM is the full width at half maximum of the first partial radial distribution function peak (estimated by a Gaussian fit) and d 0 corresponds to the mean interatomic distance (calculated as the maximum of the first peak), i.e. d 0 = 5.95 Å for Zn-Zn and d 0 = 2.0 Å for Zn-N.

The coordination number for nitrogen atoms around the zinc cation is computed by taking a cut-off radius of 2.5 Å, a value chosen from the Zn-N partial radial distribution function at room temperature. We checked that the precise value used does not influence the outcome of the calculations, nor does the choice of a discontinuous criterion (vs. the use of a damping function near the cut-off value).

In order to compute the total porous volume, we used the freely available software Zeo++. [START_REF] Pinheiro | Characterization and comparison of pore landscapes in crystalline porous materials[END_REF][START_REF] Martin | Addressing challenges of identifying geometrically diverse sets of crystalline porous materials[END_REF][START_REF] Willems | Algorithms and tools for high-throughput geometry-based analysis of crystalline porous materials[END_REF] It uses a geometric decomposition of space to compute the accessible and non-accessible volume to a sphere of a given radius. We have taken a value a 1.2 Å simulating the porous volume as seen by a helium molecule, calculating the distribution of instantaneous total pore space (sum of the accessible and the non-accessible volume) along the MD trajectories at each temperature.

Results

Structural changes upon heating and melting

In order to understand the structural evolution and onset of disorder of each framework associated with an increase of temperature, we first look at the Zn-N radial distribution functions (RDF)

for each system upon heating (Figure 1, top panel). In all cases, the overall trend observed is the same: the peaks around 2 Å and 4 Å widen, but remain present, while the other peaks flatten and ultimately disappear at very high temperature. Another common feature of interest is the fact that, above 800 K, the RDF does not go to zero between the first two peaks, indicating the presence of nitrogen exchanges between the first to the second coordination spheres of the zinc cations. Beyond this, we also observe differences between the different frameworks, at low and intermediate temperatures. There, ZIF-8 appears more structured than ZIF-zni, which better defined and narrower peaks; and ZIF-zni itself shows more structuration than ZIF-4. Moreover, the liquid obtained from ZIF-4 has features that are closer to its crystalline precursor than the two other frameworks. These trends observed on the partial Zn-N radial distribution functions are confirmed when we look at the total RDF (Figure S1), containing all pair of atoms, showing the crucial role of the zinc-nitrogen coordination bond in monitoring and understanding the process of melting in ZIFs.

To quantify this analysis, we computed the Lindemann ratio for the three ZIFs as a function of temperature (Figure 1, bottom panel). Computed from the width w of the first peak at distance d 0 in the different partial radial distribution function, the Lindemann ratio ∆ = w/d 0 is as a measure of the fluctuation of interatomic distances. A jump in the Lindemann ratio, which is rather linear at low temperature, and the crossing of the threshold value of 15%, are typically understood to evidence the melting of a solid. This occurs in the temperature range of 1,200 K to 1,500 K for ZIF-4, in the range of 1,500 K to 1,750 K for ZIF-zni, and above that value for ZIF-8. We can thus assign an "order" to the melting temperatures of the three frameworks, with T m (ZIF-4) < T m (ZIF-zni) < T m (ZIF-8). This may explain why ZIF-8 is experimentally observed not to melt, because its melting temperature is too high compared to its decomposition temperature.

Finally, we computed the average mean square displacement per atom (summed over all atoms of a framework), over a reference length of time (τ = 75 ps). We observe on Figure S2 that diffusion really starts between 1,250 K and 1,500 K for ZIF-4 and only above 1,500 K for ZIF-8. Diffusion in ZIF-zni seems to be hindered by its very dense structure, qualitatively different from the porous structures of ZIF-4 and ZIF-8. As a result, diffusion only appears above 1,750 K in this framework, whereas its Lindemann ratio at this temperature is higher than the one for ZIF-8. The comparison between the three frameworks indicates again that ZIF-8 ordering is more affected when melting.

Although the "theoretical" melting temperatures we report are well above experimental melting temperatures (around 840 K for ZIF-4 [START_REF] Bennett | Hybrid glasses from strong and fragile metal-organic framework liquids[END_REF] ), we believe our modelling still provides insights on the underlying mechanism happening at the "real" melting temperatures. Ab initio methods do not allow long time scales and very large systems, but by going to higher temperatures and checking the Arrhenian behaviour of the rare events observed, such as bond cleavages, we are able to compare between these three frameworks.

Thermodynamics of melting

In order to better understand the thermodynamics of ZIF melting, we proceeded used the Zn-N distance as a reaction coordinate, following the breaking and reformation of zinc-imidazolate coordination. From the partial radial distribution function between Zn and N, g Zn-N (r), we computed the potential of mean force (PMF), F(r) = -k B T ln g(r), at all temperatures. The resulting free energy profiles are shown in Figure 2, and present similar features: the barrier at 300 K is too high to be measured (since no dissociation event is observed during the time of our simulations), and lower with increasing temperature. The barrier height, when it can be measured, is the activation free energy ∆F ‡ associated with the breaking of the Zn-N coordination bond. Its temperature dependence is displayed as a van 't Hoff plot in Figure S3, and a linear approximation (van 't Hoff law, ∆F ‡ (T ) = ∆U ‡ -T ∆S ‡ ) is used in each case -with a reasonable fit -to obtain values for the activation energy and entropy, ∆U ‡ and ∆S ‡ respectively. The values obtained are given in Table 1.

In all cases the free energy barrier is dominated by the energetic term, while the activation entropy brings a small stabilization at high temperature. We observe that the free energy barriers encountered at the estimated melting temperatures are of 6 k B T for ZIF-8 (at 1,500 K), and 8 k B T

for ZIF-4 and ZIF-zni (at 1,200 K), significantly higher than thermal fluctuations. Thus, we confirm for ZIF-8 and ZIF-zni what was found for ZIF-4, namely that melting occurs through an activation process leading to rare events of bond breaking. [START_REF] Samanta | Microscopic mechanisms of equilibrium melting of a solid[END_REF] That is, melting occurs not at the limit of stability of the solid phase, but at a point where it is still metastable.

Comparing the potential of mean force for the Zn-Im coordinate (where Im is the center of mass of the imidazolate linker) allows us to distinguish better between the behavior of the different frameworks. In fact, as shown in Table 1, the thermodynamic parameters obtained from Zn-N to Zn-Im potential of mean forces are not the same for the three frameworks. For ZIF-4, there is almost no variation with less than 3% difference for the enthalpic and entropic terms. ZIF-zni shows a similar behaviour with less than 3% difference for both terms too. Although these two frameworks show that Zn-N bond breaking can almost be identified to Zn-Im bond breaking in terms of temperature dependence, the story for ZIF-8 is pretty different. In fact, the increase in enthalpy is 43% and even 60% for the entropic term. It explains why ZIF-8 seems to melt at a higher temperature than ZIF-4 as the resulting activation free energy with this reaction coordinate is still equivalent 11.5 k B T at 1,200 K and only goes down to 7 k B T at 1,500 K. Moreover, we observe that there is a large difference between the ∆U ‡ of the Zn-N coordinate and that of the Zn-Im coordinate in ZIF-8. This shows that the difference in behavior between ZIF-8 and the two other frameworks does not find its root in the intrinsic strength of the Zn-N bond itself (as they are chemically equivalent), but in the influence of the environment of the framework.

Microscopic mechanism

As we have seen, melting appears through creation of defects in the coordination structure, where the Zn-N coordination bonds are broken and later re-formed. In our molecular dynamics simulations, we can track directly the formation of such defects by computing the Zn coordination number as a function of framework and temperature. These averaged coordination numbers are shown in an Arrhenius plot on Figure 3. This leads to energies ε associated with defect formation, which are 71 kJ.mol -1 for ZIF-8, 56 kJ.mol -1 for ZIF-4, 67 kJ.mol -1 for ZIF-zni. These values are in the same range of magnitude, but we note here again that ZIF-4 exhibits the smallest defect formation energy, in line with thermodynamic and structural properties.

Figure 3 shows the proportion of n-fold zinc atoms (n = 0 to 5) for the three ZIFs studied, as functions of temperature. For all structures, we see that increase in temperature leads to undercoordination of the zinc atoms through breaking of the Zn-N bonds. The trend we see on the four-fold coordination matches the structural quantities studied before, as the under-coordination comes at lower temperatures for ZIF-4 (35% drop between 1,200 K and 1,500 K) but is then much larger for ZIF-8 (66% drop between 1,500 K and 1,750 K), consistently with a greater disorganization due to diffusion. ZIF-zni four-fold coordination really jumps at higher temperatures, with a 45% drop between 1,750 K and 2,000 K. The microscopic mechanism, which we first unveiled in ZIF-4, is broadly similar in ZIF-8 and ZIF-zni. The transition into the liquid state involves rapid, activated events of imidazolate linker switching from a zinc cation. From an initial fourfold coordinated zinc, one imidazolate linker moves away and is then replaced by a neighbouring Im group with a dangling nitrogen lone pair.

We analyzed these "linker exchange" events, which are depicted in Figure 4. Inspired by Laage et al.'s statistical treatment of the molecular mechanism of reorientation in liquid water [START_REF] Laage | A molecular jump mechanism of water reorientation[END_REF][START_REF] Laage | On the molecular mechanism of water reorientation[END_REF] , we averaged over all events during a molecular dynamics trajectory to plot the average distance of the outgoing and incoming N atoms during an exchange-taking the reference time t = 0 when both nitrogen atoms are at equal distance from the zinc ion. The resulting average dynamic, shown on Figure 4, show that the bond cleavage mechanism is similar for the three frameworks. The only small difference is the coordination state during the exchange for ZIF-8. Although for ZIF-zni, as for ZIF-4, the exchange happens under slight undercoordination ( -2%), in the case of ZIF-8 it actually happens under a small overcoordination ( 3%). However, this can be explained by the fact that the average coordination number in ZIF-8 liquid is indeed much lower ( 2.3) than it is in ZIF-4 and ZIF-zni ( 3.3).

Presented in Table 2 are the activation energies associated with the Zn-N and Zn-imidazolate cleavage frequencies. ZIF-zni shows the highest value (137 kJ.mol -1 ) in terms of activation energy for Zn-N bond breaking which confirms the difficulty to thermodynamically initiate the melting in this framework. In ZIF-4 the two processes were clearly related as the difference in activation energy was only 6.5% more for Zn-Im than for Zn-N. However for ZIF-8 this difference rises to 58%. This can be understood by looking at the last column in Tables S1, S2 and S3. Indeed the proportion of Zn-N cleavage corresponding to a Zn-Im/mIm cleavage is slowly growing in ZIF-4 while it goes from 18% to 94% between 1,500 K and 2,000 K in ZIF-8. ZIF-zni does show a 28% increase between the two activation energies, which can be ascribed to the imidazolate diffusion being hindered due to its denser structure.

Characterization of the liquid ZIF

We have characterized the liquids obtained at temperatures above the melting point by computing, from our FPMD simulations, the translational diffusion coefficients of zinc cations and imidazo-late groups (see Table S4). The behaviour observed is qualitatively similar to what we saw with the global RMSD presented in Figure S2: diffusion in ZIF-8 starts to be significant at higher temperature than in ZIF-4. For instance at 1,500 K, D Im is equal to 7.7 10 -10 (m 2 .s -1 ) which is more than ten times the diffusion coefficient of methylimidazolate in ZIF-8 (4. 7 10 -11 ) at the same temperature, but this ratio is inversed at 2,000 K where D ZIF-4

Im is about one third of D ZIF-8 mIm . ZIF-zni shows smaller diffusion coefficients all the way as expected. Nonetheless, the three ZIFs exhibit the same symmetry between the transport of zinc and imidazolate groups, as their diffusion coefficients are always very close, a behaviour commonly encountered in classical ionic liquids with strong association.

From the diffusion coefficients at various temperatures, we fitted Arrhenius laws -although we note that there are clearly uncertainties and fluctuations in the values of diffusion coefficients, due to the relatively short trajectories allowed on these large systems by first-principles methods.

Characteristic activation energies for all frameworks are presented in Table S4. For example, in the case of Zn we have the following values: 302 kJ.mol -1 for ZIF-8, 105 kJ.mol -1 for ZIF-4 and 284 kJ.mol -1 for ZIF-zni. The ordering of activation energies associated with diffusion match the one associated with breaking the Zn-imidazolate, thus confirming a diffusive mechanism. Contrary to fully dissociated ionic liquid though, zinc atoms and imidazolate ligands stay strongly associated in average as the aforementioned coordination numbers show.

Finally, when talking about MOFs, porosity is one of the most crucial property that comes to mind. We have thus investigated the porosity of our frameworks as a function of temperature.

For that matter, we performed statistical analyses using a probe with a 2.4 Å diameter, the kinetic diameter of helium. The results are presented in Figure 5. ZIF-zni is a dense crystalline structure, with almost no porosity, and as expected this changes very little with temperature (Figure S4). The ZIF-4 liquid, as we demonstrated in previous work, retains most of the modest porosity of its parent crystalline phase, with instantaneous "pockets" or voids in the liquid, even at high temperatures. In stark contrast, ZIF-8 -well known for high porosity of its crystal structure -sees its pore space suffer a lot from the melting process. More than half of its porosity is lost, and the distribution of instantaneous porous volume is rather broad. We note that it is, however, still larger than that of the ZIF-4 liquid.

Conclusions and perspectives

We recently introduced the general term of "MOF liquid", [START_REF] Gaillac | Liquid metal-organic frameworks[END_REF] following a combined experimental and computational study of the melting of ZIF-4, using in situ variable temperature X-ray, ex situ neutron pair distribution function experiments, and first-principles molecular dynamics simulations. This initial study paved the way for investigations into the questions it had opened, about the generality of the phenomenon and the properties of the MOF liquids. Here, we have used first-principles molecular simulations in order to investigate the question of transferability of the melting phenomenon among three different ZIF frameworks, and the influence of the nature of the parent crystalline ZIF on the physico-chemical properties of the resulting liquid. We conclude that the phenomenon is general, and also that the mechanism by which melting occurs is similar in all three frameworks. We also show that there is a clear impact of the characteristics of the ZIF crystals on both the melting, with ZIF-8 melting above its experimental temperature of decomposition, and thus not a candidate. We link this to the high free energy barrier for the detachement of imidazolate linkers from the metal cations. We hypothesize that this may be due to the high porosity of ZIF-8:

while for denser ZIFs, like ZIF-4 and ZIF-zni, the movement of the ligand allows it to keep good contact with other ligands and be stabilized by dispersive interactions, in highly porous ZIF-8 the detachement of the imidazolate group implies that it be very isolated in the intermediate state, and thus a high energetic barrier.

Thus, framework topology plays an important role in the determination of melting, as well as in the properties and retention of the porosity in the liquid phase. There is thus a need to study the melting of other ZIF materials, combining experimental and theoretical forces, in order to confirm this hypothesis and establish general criteria for low melting temperatures in ZIFs. This will require the study of ZIFs with different organic ligands, different topologies, but also different metal cations (such as cadmium) [START_REF] Baxter | A comparison of the amorphization of zeolitic imidazolate frameworks (ZIFs) and aluminosilicate zeolites by ball-milling[END_REF] that can influence the strength of the zinc-imidazolate interaction liquids with tunable chemical and physical properties, as is done for crystalline MOF architectures that can be designed and tailored for specific applications. Such materials would be of interest for liquid phase separations, homogeneous catalysis and ion transport. They could also be used as intermediaries to obtain mechanically and thermally stable porous MOF glasses.

  , and a relative cutoff (keyword REL_CUTOFF in CP2K) of 40 Ry for high accuracy. Valence electrons were described by double-zeta valence polarized basis sets and norm-conserving Goedecker-Teter-Hutter 25 pseudopotentials all adapted for PBE (DZVP-GTH-PBE) for H, C and N or optimized for solids (DZVP-MOLOPT-SR-GTH) in the case of Zn. The unit cell studied for ZIF-4 (space group Pbca) is the orthorhombic primitive unit cell, which contains 272 atoms, with cell parameters a = 15.423 Å, b = 15.404 Å, c = 18.438 Å, and α = β = γ = 90°. The unit cell studied for ZIF-8 (space group I -43m) is the orthorhombic primitive unit cell, which contains 276 atoms, with cell parameters a = b = c = 16.993 Å, and α = β = γ = 90°. The unit cell studied for ZIF-zni (space group P1) is the primitive unit cell, which contains 272 atoms, with cell parameters a = 17.570 Å, b = 17.613 Å, c = 17.613 Å, and α = 138.433°, β = 97.256°, γ = 97.257°. Representative input files for the molecular dynamics simulations are available as supporting information, and online in our data repository at https: //github.com/fxcoudert/citable-data.

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Top panel: Radial distribution functions for the pair of zinc and nitrogen atoms in the three ZIFs as a function of temperature. Bottom panel: Comparison of the evolution of the Lindemann ratios for the 3 ZIFs showing melting is happening at different temperatures.

Figure 4 .

 4 Figure 4. Top panel: Behaviour during an exchange of a nitrogen atom by another nitrogen atom, in the first coordination sphere of a zinc cation, averaged over all such events (all exchanges on all zinc cations) in the three ZIFs. The distance between the incoming nitrogen and the zinc is plotted in red, and that between the outgoing nitrogen and the zinc in blue. The green curve corresponds to the average coordination number of the zinc cation involved in the exchange. The flat dashed lines are the average coordination number over the whole simulation at 300 K (light blue) and 2,000 K (green). Bottom panel: Visualization of a representative imidazolate exchange event. Zn: green, N (initially coordinated): light blue, N: blue, N (coordinated after exchange): orange, C: grey. Bottom panel reproduced from Ref. 15, copyright 2017, courtesy of François-Xavier Coudert.

Figure 5 .

 5 Figure 5. Comparison of temperature-dependent porosity between ZIF-4 (green) and ZIF-8 (red).The relative stability of the porous volume in ZIF-4 is striking compared to the collapse of the one in ZIF-8.

Table 1 .

 1 The first column shows the activation enthalpy associated with a bond cleavage, the second column shows the activation entropy and the third column shows the extrapolated activation free energy at the experimental melting temperature for ZIF-4 of 840 K.

	8.

Table 2 .

 2 This table presents the activation energies needed to break a Zn-N bonds (first column)and Zn-Im/mIm bonds (second column) in the three ZIFs. volumes, crystalline structures. Additional statistics on the mechanism of bond cleavage and diffusion coefficients.

Acknowledgements

We thank Anne Boutin, Alain Fuchs, and Thomas Bennett for fruitful discussions. This work benefitted from the financial support of ANRT (thèse CIFRE 2015/0268) and access to HPC platforms provided by a GENCI grant (A0030807069).

Supporting Information

This material is available free of charge via the Internet at http://pubs.acs.org.

Total radial distribution function, mean square displacements, activation free energies, porous


Laura Scalfi

Guillaume Fraux

Anne Boutin
email: anne.boutin@ens.fr

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

†pasteur

Structure and Dynamics of Water Confined in Imogolite Nanotubes

We have studied the properties of water adsorbed inside nanotubes of hydrophilic imogolite, an aluminum silicate clay mineral, by means of molecular simulations. We used a classical force field to describe the water and the flexible imogolite nanotube, and validated it against data obtained from first-principles molecular dynamics. With it, we observe a strong structuration of the water confined in the nanotube, with specific adsorption sites and a distribution of hydrogen bond patterns. The combination of number of adsorption sites, their geometry and the preferential tetrahedral hydrogen bonding pattern of water leads to frustration and disorder. We further characterize the dynamics of the water, as well as the hydrogen bonds formed between water molecules and the nanotube, which are found to be more than one order of magnitude longer than water-water hydrogen bonds.

1

Introduction

Imogolite is an aluminum silicate clay mineral, with formula Al 2 SiO 3 (OH) 4 , that was first discovered in volcanic ashes in Japan. [START_REF] Yoshinaga | Imogolite in some ando soils[END_REF] It is the only known aluminosilicate material that spontaneously forms inorganic nanotubes, and there exists no planar equivalent material in nature -unlike carbon nanotubes, of which graphene is the planar form. Imogolite nanotubes are monodisperse in diameter, and can be readily synthesized with controlled length and diameter, for example by substituting silicon with germanium 2 . This precise control of the imogolite nanotube dimensions is interesting for applications that rely on onedimensional pore channels, in fields such as nanofluidic devices, membranes for filtration and separation, desalination, etc. Moreover, from a theoretical point of view, its hollow cylindrical topology and tuneable size make it a very attractive model to study the properties of fluids under confinement.

Although imogolites have been the subject of several experimental studies, only a few theoretical works have been carried out until now. In particular, there is relatively little data -experimental or computational -on the hydration of these nanotubes and the behavior of the confined water molecules inside the pore space. Imogolite nanotubes differ markedly from the more common carbon nanotubes in both geometry and chemistry, and the behavior of water inside their pores is still very much open. It is experimentally challenging to address these issues, especially differentiating the water inside the nanotubes from the water outside, and a computational approach through molecular modeling is therefore a natural complement to the published experimental results.

A single imogolite nanotube presents two surfaces, both available for adsorption. The outer surface is composed of a gibbsite-like sheet of aluminum octahedra. The inner surface is formed of silicon tetrahedra with one hydroxyl group pending and the three other corners of the tetrahedron linked to three aluminum octahedra. Therefore both surfaces are covered with hydroxyl groups that are expected to have a hydrophilic behavior. Based on electrostatic calculations, authors have suggested that the inner surface is strongly hydrophilic whereas the outer surface might be more hydrophobic. [START_REF] Gustafsson | The Surface Chemistry of Imogolite[END_REF][START_REF] Guimarães | Imogolite Nanotubes: Stability, Electronic, and Mechanical Properties[END_REF] Most of the theoretical studies on imogolite nanotubes focus on the energetics of empty nanotubes. The sharp monodispersity in nanotube diameters has been explained by both quantum chemistry and classical studies, that show the strain energy has a minimum for a given nanotube diameter, contrary to carbon nanotubes where the energetically favorable structure is that of infinite diameter (i.e., the graphene slab [START_REF] Guimarães | Imogolite Nanotubes: Stability, Electronic, and Mechanical Properties[END_REF][START_REF] Tamura | Molecular Dynamics Modeling of Tubular Aluminum Silicate: Imogolite[END_REF][START_REF] Zhao | Energetic Minimum Structures of Imogolite Nanotubes: A First-Principles Prediction[END_REF][START_REF] Demichelis | Structure and energetics of imogolite: a quantum mechanical ab initio study with B3LYP hybrid functional[END_REF] ). Several studies also focused on computing vibrational spectra [START_REF] Tamura | Molecular Dynamics Modeling of Tubular Aluminum Silicate: Imogolite[END_REF][START_REF] Alvarez-Ramírez | Ab initio simulation of the structural and electronic properties of aluminosilicate and aluminogermanate nanotubes with imogolite-like structure[END_REF][START_REF] Konduri | Controlling Nanotube Dimensions: Correlation between Composition, Diameter, and Internal Energy of Single-Walled Mixed Oxide Nanotubes[END_REF] and studying the energetics and dynamics of the rolling of the nanotubes themselves. [START_REF] Lee | Origin of the Strain Energy Minimum in Imogolite Nanotubes[END_REF][START_REF] González | Model for Self-Rolling of an Aluminosilicate Sheet into a Single-Walled Imogolite Nanotube[END_REF][START_REF] González | Mechanical Response of Aluminosilicate Nanotubes under Compression[END_REF] Finally, other works have pointed at more complex aspects, such as defects [START_REF] Gustafsson | The Surface Chemistry of Imogolite[END_REF] and deformation or hexagonalisation of the nanotubes. [START_REF] Tamura | Molecular Dynamics Modeling of Tubular Aluminum Silicate: Imogolite[END_REF][START_REF] Amara | Hexagonalization of Aluminogermanate Imogolite Nanotubes Organized into Closed-Packed Bundles[END_REF][START_REF] Creton | Molecular dynamics study of hydrated imogolite : 2. Structure and dynamics of confined water[END_REF] Because adsorption on the outer surface of the nanotubes depends strongly on experimental conditions affecting the bundling of the nanotubes, and therefore the spacing between them, we will focus in the present work on the adsorption and behavior of water confined inside imogolite nanotubes. Moreover, the strong curvature and limited space inside the nanotube provide for strong confinement effects on which little experimental data is available.

We have used a classical force field to describe the water and the flexible imogolite nanotube, which we validate against structural data obtained from first-principles molecular dynamics.

From classical grand canonical Monte Carlo and molecular dynamics simulations, we observe a strong structuration of the water confined in the nanotube, with specific adsorption sites and a statistical distribution of hydrogen bond patterns. The combination of number of adsorption sites, their geometry and the preferential tetrahedral hydrogen bonding pattern of water leads to frustration and disorder. We also characterize the dynamics of the water, as well as the hydrogen bonds formed between water molecules and the nanotube.

Systems

Imogolite structure

Imogolites are aluminum silicate clay minerals with chemical formula (OH) 3 Al 2 O 3 SiOH.

Their structure was first described by Cradwick et al. [START_REF] Cradwick | a Hydrated Aluminium Silicate of Tubular Structure[END_REF] from electron diffraction measurements as a cylindrical assembly of silicon tetrahedra and aluminum octahedra (see Figure 1). This sheet spontaneously folds into a nanotube because of bond length mismatch and formation of hydrogen bonds. [START_REF] Lee | Origin of the Strain Energy Minimum in Imogolite Nanotubes[END_REF][START_REF] González | Model for Self-Rolling of an Aluminosilicate Sheet into a Single-Walled Imogolite Nanotube[END_REF] The imogolite nanotubes are characterized using the same nomenclature as carbon nanotubes. Both natural and synthetic imogolite nanotubes exhibit a zig-zag folding, with a variable number 𝑛 of (OH) 3 Al 2 O 3 SiOH units (called gibbsite units) along the nanotube circumference. The cylindrical unit cell then contains 2𝑛 gibbsite units.

Natural imogolite is a nanotube with size 𝑛 = 12, while depending on the details of the synthesis conditions, synthetic imogolites correspond to values of 𝑛 typically between 12 and 14, and sometimes larger.

In order to produce imogolite nanotube models, we started from the structure of Cradwick et al. [START_REF] Cradwick | a Hydrated Aluminium Silicate of Tubular Structure[END_REF] for a 𝑛 = 10 nanotube with 𝐶 2𝑛 symmetry. We formed a flat gibbsite-like sheet by unfolding the nanotube and adding hydrogen atoms (not present on the electron diffraction data). We relaxed the slab geometry and unit cell with density functional theory (DFT, see computational details in next section) in the 𝐶 𝑚 space group. From this relaxed planar structure, we rolled back a nanotube of size 𝑛 = 12. Using the Bilbao Crystallographic Server [START_REF] Aroyo | Bilbao Crystallographic Server: I. Databases and crystallographic computing programs[END_REF] , we identified the space groups of the nanotube to be 𝑃 4/𝑚. We relaxed the nanotube structure and cell parameters with periodic DFT calculations, with large inter-nanotube spacing (60 Å) so that there are no interactions between nanotubes. We show on Figure 1, two views from top and side of the 𝑛 = 12 nanotube with three unit cells along the 𝑧 axis.

The internal surface of an imogolite nanotube can be described as a periodic sequence of silicon rings along the nanotube axis. Two adjacent rings are rotated of 𝜋/𝑛. When dry, silanol groups within a ring form hydrogen bonds so that all the hydroxyl groups are in a plane normal to the nanotube axis, as shown in Figure 1. The internal diameter, computed between internal oxygens (noted O int ), is 12.8 Å and the external diameter, computed between external oxygens (noted O ext ), is approximately 22 Å. Konduri et al. [START_REF] Konduri | Water in Single-Walled Aluminosilicate Nanotubes: Diffusion and Adsorption Properties[END_REF] highlighted that as a consequence the inner pore is not smooth and uniform as for the carbon nanotube, but the pore is rugged: narrower where there are silanol groups and larger between these silicon rings. One can picture the internal cavity as a hollow cylinder with equally spaced circular furrows in the circumference.

Nanotubes packing

Imogolite nanotubes can have a length from 10 nm up to a few micrometers and experiments show that they tend to pack into bundles when dry. The periodic pattern seems to be a monoclinic assembly [START_REF] Ackerman | Gas/vapor adsorption in imogolite: a microporous tubular aluminosilicate[END_REF][START_REF] Mukherjee | Phenomenology of the Growth of Single-Walled Aluminosilicate and Aluminogermanate Nanotubes of Precise Dimensions[END_REF] with 𝛽 ≈ 78°. Given our focus on adsorption inside the nanotubes, following other theoretical studies, we chose to represent the nanotubes in a hexagonal packing that is close to the monoclinic one. [START_REF] Alvarez-Ramírez | Ab initio simulation of the structural and electronic properties of aluminosilicate and aluminogermanate nanotubes with imogolite-like structure[END_REF][START_REF] Konduri | Water in Single-Walled Aluminosilicate Nanotubes: Diffusion and Adsorption Properties[END_REF][START_REF] Zang | Flexibility of Ordered Surface Hydroxyls Influences the Adsorption of Molecules in Single-Walled Aluminosilicate Nanotubes[END_REF] Lattice constants are chosen to ensure close contact (but no overlap) between neighboring nanotubes, at 𝑎 = 𝑏 = 24.2 Å. This choice is consistent with experiments and other theoretical studies. Along the tube axis, we studied supercells with lattice parameters 𝑐 = 𝑛 × 8.486 Å where 𝑛 = 1, 3 or 5 (where the value of 8.486 Å was determined by geometry optimization with variable cell).

Simulation methods

We performed simulations at two different levels of theory: at the quantum chemistry level, we performed Density Functional Theory calculations for both energy minimization as well as first principles molecular dynamics. We also used classical force field-based simulations, for both classical molecular dynamics and Monte Carlo simulations.

Density Functional Theory calculations

Energy minimization of the nanotubes and slabs were performed with Density Functional Theory (DFT) calculations with periodic boundary conditions and the use of symmetry to reduce computational cost. We used these calculations with the CRYSTAL14 21 software package. The exchange-correlation functional used was the solid-adapted generalized gradient approximation (GGA) PBEsol, and the basis set for all atoms (Si, Al, O, H) was double-𝜁 valence polarized (DZVP).

We also performed DFT-based molecular dynamics simulations, aka first principles or ab initio molecular dynamics (MD), as a way to generate reference data for comparison and validation of classical force fields. These simulations were run with CP2K. [START_REF] Hutter | 2K: atomistic simulations of condensed matter systems[END_REF] Temperature was kept constant at 300 K with a CSVR (canonical sampling through velocity rescaling) thermostat, [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] with a time constant of 200 fs. The timestep was 0.5 fs, and the system was deuterated for computational convenience. After 5 ps of equilibration, production trajectories were accumulated for 30 ps.

Classical simulations

Although DFT is a good model to extract the structure and short-time dynamics of the system, simulations longer than a few tens of picoseconds are not readily accessible. This is why we used classical simulations, based on a force field approximation of the intra-and inter-molecular interactions, to study structure and dynamics of water in imogolites at larger timescales. We ran MD simulations in the canonical (𝑁, 𝑉, 𝑇 ) ensemble with the LAMMPS code, 24 using a timestep of 0.5 fs. After 100 ps of equilibration, we collected trajectories from 200 ps to 50 ns, depending on the properties studied in each case.

We also performed series of Grand Canonical Monte Carlo (GCMC) simulations with fixed chemical potential, volume and temperature (𝜇, 𝑉, 𝑇 ), to simulate the behavior of the system on equilibrium with a reservoir at a given pressure linked to the chemical potential.

In the gas phase, water was considered to be an ideal gas, and the chemical potential is then easily linked to the gas pressure. [START_REF] Desbiens | Water Condensation in Hydrophobic Silicalite-1 Zeolite: A Molecular Simulation Study[END_REF] Series of GCMC simulations were used to compute adsorption isotherms, with water vapor pressure going from 2 to 3 600 Pa.

To describe the interactions of the nanotube, we relied on the CLAYFF force field, [START_REF] Cygan | Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field[END_REF] which has been extensively used in the literature. [START_REF] González | Mechanical Response of Aluminosilicate Nanotubes under Compression[END_REF][START_REF] Konduri | Water in Single-Walled Aluminosilicate Nanotubes: Diffusion and Adsorption Properties[END_REF][START_REF] Zang | Flexibility of Ordered Surface Hydroxyls Influences the Adsorption of Molecules in Single-Walled Aluminosilicate Nanotubes[END_REF] Water was described in the flexible Single-Point Charge (fSPC) model, [START_REF] Teleman | A molecular dynamics simulation of a water model with intramolecular degrees of freedom[END_REF] which is naturally suited for coupling with CLAYFF.

CLAYFF is a general force field developed to model clay minerals, that have the same chemical nature of the imogolite. It relies almost exclusively on non-bonded interactions, where interatomic potentials are the sum of electrostatic (Coulombic) and dispersive-repulsive (Lennard-Jones 12-6) interactions. In addition to these non-bonding interactions, CLAYFF includes a harmonic bond term for hydroxyl bonds (O-H stretching). Moreover, CLAYFF can be extended by the use of an additional M-O-H harmonic bending potential, although this term is not generally used and not necessary to reproduce clay structures. We call the resulting force field "extended CLAYFF", and in this work we have used both CLAYFF and extended CLAYFF. We show later that the M-O-H bending potential is crucial for simulating imogolite nanotubes, especially for adsorption properties.

Results and discussion

Validation of the force field

The CLAYFF force field has been used in the literature for classical molecular simulations of imogolite nanotubes, [START_REF] González | Mechanical Response of Aluminosilicate Nanotubes under Compression[END_REF][START_REF] Konduri | Water in Single-Walled Aluminosilicate Nanotubes: Diffusion and Adsorption Properties[END_REF][START_REF] Zang | Flexibility of Ordered Surface Hydroxyls Influences the Adsorption of Molecules in Single-Walled Aluminosilicate Nanotubes[END_REF] owing to the similarities between clays and imogolite. Here, we set out to check its validity in the description of both neat imogolite nanotubes as well as water-filled tubes. As a reference, we compared it to trajectories obtained from first-principles molecular dynamics, where the evaluation of interatomic forces is done at the DFT level. We found that the structure of the nanotube itself is relatively rigid, and dictated by equilibrium In contrast, augmenting the CLAYFF force field with the additional M-O-H bending term (the "extended CLAYFF") leads to a clear improvement of the imogolite description.

The distributions of both SiOH and AlOH angles are closer to those obtained from firstprinciples MD, although the match is not perfect -as could be expected for a generic force field without ad hoc reparameterization. The average value of the SiOH angle with extended CLAYFF is 122°, which is 7°larger than the reference data. The same is true of the outer AlOH groups.

These discrepancies in imogolite structure strongly affect water-imogolite interactions, and from here water adsorption, structure and dynamics. To further assess the necessity of using the extended CLAYFF instead of CLAYFF, we performed simulations of fully hydrated nanotubes, and compared them against data from first-principles MD. On the right panel of Figure 2, we show the density profiles along the nanotube axis 𝑧 for the internal hydrogens.

First, we observe that in the empty nanotube, the density distributions are all around 𝑧 = 0 (the plane in which are located the silicon atoms). This allows hydrogen-bonding with other silanol groups in the same 𝑧 = 0 plane. Secondly, upon hydration, there is a clear move of the SiOH hydrogen atoms toward symmetric positions at ±0.75 Å, due to hydrogen bonding with the water molecules -a detailed analysis of this hydrated structure will be given below. This change in the H localization, which we observe in the first-principles MD, is well reproduced by the extended CLAYFF force field. However, it is markedly different in regular CLAYFF, where we see a unique broad population with two bumps. Thus the distribution of internal hydrogen atoms, fundamental for adsorption studies and the description of hydrated imogolites, is best described by the extended CLAYFF model. This aspect has been recently studied by Pouvreau et al. [START_REF] Pouvreau | Structure of Hydrated Gibbsite and Brucite Edge Surfaces: DFT Results and Further Development of the ClayFF Classical Force Field with Metal-O-H Angle Bending Terms[END_REF] that similarly showed the necessity of a M-O-H bending term to properly describe basal and edge surfaces of gibbsite and brucite. From here on, we will use exclusively the extended CLAYFF force field.

Water adsorption

In order to characterize water adsorption in the gas phase and obtain physically meaningful water uptake values, we performed GCMC simulations for water pressure between 2 and 3600 Pa. Full details of the simulations, done with two different water models (SPC and TIP4P05) are presented in Supporting Information, as well as isotherms. For computational reasons, rigid water models were preferred in the GCMC calculations. We observe (Figure S1) that the imogolite nanotube is hydrophilic, i.e. the filling occurs at a pressure below the bulk saturation pressure. This is in agreement with the affinity of water for silanol-rich surfaces, or zeolites with many silanol defects. [START_REF] Coudert | Water nanodroplets confined in zeolite pores[END_REF] The value of the pore filling pressure varies strongly with the nature of the model chosen for the description of the water, and the movement allowed for the silanol groups -consistent with the findings of Zang et al. [START_REF] Zang | Flexibility of Ordered Surface Hydroxyls Influences the Adsorption of Molecules in Single-Walled Aluminosilicate Nanotubes[END_REF] . We find transition pressures of 0.1 kPa for TIP4P05, and 1 kPa for the SPC water model. The order of magnitude of those values are within the range of experimental data, for example the gravimetric study of Konduri et al. [START_REF] Konduri | Water in Single-Walled Aluminosilicate Nanotubes: Diffusion and Adsorption Properties[END_REF] , which features a smooth isotherm with water uptake in the same range (0.1 to 1 kPa; see Figure S1). In contrast to the adsorption pressure model dependence, there is a excellent agreement between the various models on the saturation uptake of water once the pore is filled, i.e. on the density of the water adsorbed inside the imogolite nanotube. Uptake is found to be 10 % by weight after the transition, both for the SPC and TIP4P05 water models.

We thus derived from the GCMC calculations an initial configuration for the simulation of the fully filled imogolite nanotube, from the plateau of the SPC adsorption isotherm. This initial configuration contains three nanotube (𝑛 = 12) unit cells along the 𝑧 axis, 98 water molecules inside the nanotubes and 18 water molecules outside. It is represented in Figure 3.

We used this fully hydrated tube as a starting configuration for classical molecular dynamics simulations, with the extended CLAYFF force field and fSPC water model.

Structure of confined water

From the analysis of the MD trajectories of the fully hydrated imogolite, we computed density profiles of all atom types (water O w and H w ; silanol groups Si, O int , and H int ). They are plotted in Figure 4 in the 𝑥𝑦 plane, i.e. as if viewed "from the top" of the (infinite) nanotube. From the water oxygen distribution O w , it is clear that there are two different populations of water molecules: strongly structured water adsorbed next to the inner surface (a first hydration layer, at distances between 4.5 and 5.8 Å from the nanotube's center), and more disordered water filling up the center of the nanotube. This is also clearly visible in the cylindrical distribution functions (Figure S2), where we see that the two populations are not completely separate, i.e. the density does not fall to zero in-between. Furthermore, we see a slight deviation of the nanotube from a purely circular form: the position of the Si atoms show a slightly hexagonal deformation, due to the symmetry of the packing of nanotubes in bundles. [START_REF] Amara | Hexagonalization of Aluminogermanate Imogolite Nanotubes Organized into Closed-Packed Bundles[END_REF] This deformation is only very small, with Si displacements of 0.1 to 0.2 Å at most, due to the rigid nature of the imogolite nanotube, especially with such a small diameter. This deformation is both linked to the hexagonal packing and the adsorption stress exerted by the water molecules -as shown previsouly for water/quartz interfaces, [START_REF] Gor | Adsorption-Induced Surface Stresses of the Water/Quartz Interface: Ab Initio Molecular Dynamics Study[END_REF] and more generally for adsorbates in "soft" porous materials. [START_REF] Neimark | Stress-Based Model for the Breathing of Metal-Organic Frameworks[END_REF][START_REF] Mouhat | Softening upon Adsorption in Microporous Materials: A Counterintuitive Mechanical Response[END_REF] In the case of imogolite, the stress exerted is counterbalanced by the relative stiffness of the nanotube, although the extent and details of the deformation depend on the packing of the nanotubes. For example, Creton et al. [START_REF] Creton | Molecular dynamics study of hydrated imogolite : 2. Structure and dynamics of confined water[END_REF] used a larger distance between the nanotubes and observed a transition between an ellipsoid shape to a more cylindrical shape when increasing water densities.

To better understand the structuration of the first adsorbed layer, we plot in Figure 5 the same densities, but this time as a function of cylindrical coordinates 𝑧 and 𝜌, where 𝜌 is the distance to the central axis of the nanotube. There, the strongly preferential location for the confined water oxygen atoms becomes clear: the water molecules sit in the wider furrows, situated in-between two rings of silanol groups (SiOH rings are situated every 4.25 Å along the 𝑧 axis, and are depicted on Figure 5 by dotted lines). This preferential localization of the adsorbed water contrasts sharply with that of water inside carbon nanotubes, whose internal surface is smoother. It was, however, demonstrated in other nanopores of small dimensions with "dangling" hydroxyl groups accessible to the water for hydrogen-bonding. [START_REF] Haigis | Investigation of structure and dynamics of the hydrated metal-organic framework MIL-53(Cr) using first-principles molecular dynamics[END_REF] These molecules are strongly hydrogen-bonded to the silanol groups, as indicated by the well-defined positions for H int atoms, which are rearranged compared to their relaxed position. Indeed, in the dry imogolite nanotube, the most stable position for the internal hydrogens is in the silanol ring plane (𝑧 = 0), where they are pointing toward a neighboring O int atom. This conformation allows the formation a relatively weak hydrogen bond, because the O int -H int • • • O int distance is large. In the hydrated tube, the hydrogen atoms' most common position is shifted out of that silanol ring plane (see Figure 5). This is reflected by a change in the O br -Si-O int -H int dihedral angle (see Figure S3), which goes from 0°to 60°;

(Å) (Å) (Å) (Å)
where O br is the oxygen atom bridging the aluminum to the silicon.

We thus propose the following model to better visualize and understand the structure of the adsorbed water. We model the inner surface of the nanotube as a periodic sequence of triangles, where each silanol group SiOH is a vertex. The Si-Si and O-O radial distribution functions show two first neighbor peaks at 4.2 and 4.7 Å for silicon, and at 3.3 and 4.6 Å for internal oxygens. This means that the triangles are isosceles with two angles of 66.5°and one angle of 47°. This pattern is depicted on Figure 6. Above the center of each triangle is a potential water adsorption site. However, analysis of the sites show that no neighboring sites in the same 𝑥𝑦 plane can be occupied at the same time, due to short-distance intermolecular repulsion of water molecules, so that at most half of the adsorption sites are occupied in the filled tube. Then, each of the hydrogen atoms of the SiOH groups that form the vertices of a triangle will point toward one of the three occupied neighboring adsorption sites.

A final visualization can be obtained by a density plot of O and H atoms, where the density is mapped against 𝑧 and a circular coordinate, as if a slice of the nanotube had been "cut and unrolled". These densities are presented in Figure 7. Red dotted lines show the SiOH rings, while orange dotted lines show the triangular mesh described above. The adsorption sites on top of each triangle are clearly visible. We further see that there is a strong anisotropy of the system, where there is possible circular movement of the water molecules in-between sites in the same 𝑧 plane, while there is no observed density crossing the silanol rings, as already observed by Creton et al. [START_REF] Creton | Molecular dynamics study of hydrated imogolite : 2. Structure and dynamics of confined water[END_REF] . This shows that there is surface diffusion of water molecules, but only in the 𝑥𝑦 plane. 1. The change in silanol-silanol hydrogen bonding confirms the structural change upon adsorption: when the nanotube is empty, approximately half of the silanol groups are linked by a hydrogen bond within silicon rings (as detected per our rather strict criteria).

When hydrated, only 2% form a hydrogen bond with another silanol group. Regarding the hydrogen bonds between a water molecule and the surface, 94% of the silanol groups donate a hydrogen bond to a water molecule and each silanol group receives on average 0.96 hydrogen bonds from water. Since a silanol group can donate only 1 hydrogen bond and can receive 1 or 2 hydrogen bonds, the bonding possibilities of the surface are extensively used. Within the nanotube, on average 57% of water molecules receive or donate at least one hydrogen bond from or to a silanol group, i.e. are in the first adsorbed layer.

To further characterize the local hydrogen bond network, we considered all 4 singular hydrogen bonding behaviors possible for a single water molecule: donating to a silanol or a water molecule, or accepting from either. Each of these singular behaviors can appear zero, 2.

once, or twice. The patterns that represent more than 5% of the occurences found over the whole trajectory are described in Table 2.

The six patterns highlighted by this analysis can be divided in two groups: (i) water molecules that do not form any hydrogen bond with a silanol group (patterns 1 and 4), that can have either 4 hydrogen bonds (tetragonal arrangement) or 3 H bonds, similarly to bulk water; and (ii) molecules whose adsorption on the inner surface of the nanotube involves at least one hydrogen bond (types 2, 3, 5 and 6). It is worth noting that the mean number of hydrogen bonds for a water molecule not bonded to the surface is 3.65, slightly higher than the bulk value for fSPC water (3.60 hydroben bonds per water molecule). The most common patterns are represented in Figure 9. For adsorbed molecules, the most common pattern is an "upright" water molecule (pattern 2, Fig. 9b), in a plane bisecting the triangular adsorption site. This water molecule receives two H bonds from two silanol groups and donates one hydrogen bond to the third silanol. A water hydroxyl group is left pointing toward the center of the nanotube, free to donate a hydrogen bond to a water molecule of the second layer.

This creates a structuration of the adsorbed water beyond the first layer, which is clearly visible in the density maps (Figures 4 and5), with a favorable position for water oxygen atoms across that pending hydroxyl group. It also defines an equilibrium distance between the first and second layer of adsorption, which corresponds to a water-water hydrogen bond length.

The second most common pattern for the first layer of water is a water molecule lying "flat" on the inner imogolite surface, donating two H bonds to neighboring silanol groups and receiving one from the third silanol neighbor. This allows the formation of a fourth hydrogen bond from a water molecule outside of the first layer. In this case the direction of the hydrogen bond is less constrained, explaining the "diffuse" part of the oxygen density in the second water layer, seen in Figure 5.

The large possibilities of organization of the first layer of adsorbed water molecules, with orientational disorder, is reminiscent of the theory of Bernal and Fowler 35 regarding ice rules. From the energetics of each individual configurations, as well as the frustration arising from constraints (no more than two H bonds can be donated or received, and no two neighboring adsorption sites can be populated at the same time), we believe that a model of the orientational disorder in that strongly bound layer of adsorbed water could be constructed, possibly explaining why the adsorption sites are only occupied at nearly 38% at saturation uptake.

Water dynamics

From the structural analysis of the water density and the hydrogen bonds in the adsorbed water, there are clear hints at the dynamic nature of the system and the lability of the water molecules. We have seen that hydrogen bond patterns are statistically distributed among several configurations, and that there is anisotropic surface diffusion in the first layer of adsorbed water (see Figure 7 and accompanying text). Analysis of trajectories shows that this surface diffusion occurs without breaking all the hydrogen bonds involved. By breaking a single hydrogen bond, an adsorbed water molecule can flip around the axis formed by the remaining silanol groups, and thereby diffuses to the adjacent site, then reforms a hydrogen bond with another SiOH group. Since the triangular sites are isosceles and the surface is curved, the flip is easier and requires less configurational change within the furrows, than crossing the silicon rings. This explains the preferential circular diffusion as opposed to To do so, we calculated autocorrelation functions of the presence of water molecules, both within the first adsorption layer and in the less-ordered "core" of the tube. Due to the relatively slow dynamics of the strongly-confined water, these quantities were calculated from 50 ns classical trajectories. For all three possible donor-acceptor hydrogen bond types, we computed autocorrelation functions which are plotted in Figure 10a. From these autocorrelation functions, the lifetime of each H bond type was extracted, and are presented in Table 3. The autocorrelation decay is not stricly exponential, because the hydrogen bonds form multiple populations (first and second layer, etc.) with different decays.

First, looking at the water-water hydrogen bonds, we see that their lifetime is roughly 2 times higher in the imogolite nanotube than in bulk water (calculated at 2.67 ps from a simulation of pure fSPC water). This slow down of the water dynamics under confinement has been well studied and explained in previous work, and is due in part to the stronger organization of the water, and in part to excluded volume effects (a water molecule finding fewer partners to "switch" a hydrogen bond to). [START_REF] Fogarty | Reorientational Dynamics of Water Confined in Zeolites[END_REF]37 Secondly, we see that hydrogen bonds between water and the silanol-rich surface have a lifetime more than an order of magnitude greater than water-water H bonds. A study on hydrogen bonding of water on the 𝛼-quartz interface -whose (001) surface is chemically similar to imogolite, in that it presents geminal silanol groups -calculated persistence times of SiOH-water hydrogen bonds between 160 and 170 fs. [START_REF] Ozkanlar | Water Organization and Dynamics on Mineral Surfaces Interrogated by Graph Theoretical Analyses of Intermolecular Chemical Networks[END_REF] The extremely large difference, and the very strong bonding of water on the imogolite inner surface, can be explained by the cooperative effect of the formation of three hydrogen bonds, on top of which is present the generic confinement slow-down.

We also investigated the residence time of water molecules in the first adsorbed layer, inside the imogolite nanotube. There are two ways to define the "first layer", by using geometrical radius, or by considering molecules that form at least one hydrogen bond with the surface. We chose the latter option, which is more chemically robust. The time autocorrelation function deriving from it is plotted in Figure 10b. It shows an exponential decay, and reaches a plateau at 57%, i.e. the average percentage of water molecules that belong to that first layer. The decay time for the autocorrelation function is found to be 410 ps. This is two or four times larger than the indicative lifetimes of hydrogen bonds with the surface (donated or received by the silanol group, respectively). This confirms the existence of an "inner surface diffusion" mechanism that allows diffusion between adsorption sites without exiting the layer itself. We can also compare this lifetime to residence times in protein adsorption sites.

For example, for myoglobin [START_REF] Makarov | Residence Times of Water Molecules in the Hydration Sites of Myoglobin[END_REF] , residence times span from 10 to 450 ps where the unusually high residence times (> 80 ps) correspond to water trapped either in the cavities inside the protein or in the grooves and concave regions.

To complete this study, we extracted diffusion coefficients from mean square displacement (Figure S4) of water molecules in the nanotube. We find a diffusion coefficient in the 𝑧 axis of 3.5 10 -6 cm 2 /s, which is three times smaller than the coefficient calculated by Zang et al. [START_REF] Zang | Self-Diffusion of Water and Simple Alcohols in Single-Walled Aluminosilicate Nanotubes[END_REF] for a similar but lower water density. Moreover, we do not see a long-time diffusive behavior in the 𝑥𝑦 plane because of the confinement.

Conclusions and perspectives

This study of the properties of water adsorbed inside hydrophilic imogolite nanotubes was motivated by the understanding of the specific behavior of water due to its tight confinement, with numerous silanol groups on the inner surface of the nanotubes. Because the strong interactions between water and imogolite, and the slow dynamics of the confined water, we relied on a classical force field in order to perform molecular dynamics studies. We validated the CLAYFF force field, and found that it had to be employed in its extended versionwith a Si-O-H bending term -in order to correctly reproduce the reference data we had produced from first-principles molecular dynamics. This is in contrast with conventional clay materials, where this term is not as important, and is ascribed to the strong curvature of the imogolite surface.

We then developed a model for the structure of adsorbed water, based on the calculated atomic densities and statistical analysis of the hydrogen bonding patterns observed between water molecules and silanol groups. We show that adsorption in the first layer vicinal to the imogolite tube can be described by an anisotropic triangular lattice of adsorption sites, which are not fully populated. We showed that interactions between the water molecules adsorbed at neighboring sites are complex, based on (i) exclusion of some nearest neighbor pairs (where the sites are too close to be occupied at the same time), and (ii) the number of hydrogen bonds donated (and accepted) by the water molecules to (and from) neighboring silanol groups. Maximizing the number of hydrogen bonds created, while adhering to these rules, creates frustration in the system and leads to the emergence of a heavily disordered state. This situation is very similar to the Bernal-Fowler rules describing the orientation of water molecules in ice. [START_REF] Bernal | A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions[END_REF] Finally, we also characterized the dynamics of the confined water molecules, as well as the hydrogen bonds. We find that in addition to the generic effect of slow-down of confined water, the strong interactions of water molecules with the silanol groups giving rise to lifetimes that compare to typical residence times of water in protein adsorption sites.

In the future, it would be interesting to extend this study of the structure of adsorbed water at varying temperature, or with changes in the imogolite nanotube size (from the value of 𝑛 = 12 repeating units used in this study). Another noteworthy perspective would be to consider the possibility of partial or total functionalization of the silanol groupsreplacing for example some of the hydrophilic -SiOH groups with hydrophobic -SiOCH 3 .

This would affect the overall affinity of the internal surface of the imogolite nanotubes for water, tuning them from very hydrophilic to potentially hydrophobic, [START_REF] Bonelli | Developments in Clay Science[END_REF] and changing the thermodynamics of adsorption. It would also introduce more frustration into the adsorbed water and change the "ice rules", affecting the structure of the confined phase. This could be studied by atomistic simulations, as was done here, or at a more fundamental scale by an Ising-type model of the adsorbed water molecules, dictated by the energies of the various possible configurations of neighbording water molecules and silanol groups.

Figure 1 :

 1 Figure 1: Imogolite nanotube (𝑛 = 12) with three unit cells along the tube axis 𝑧. Color code: Aluminum (pink), Silicon (yellow), Oxygen (red), Hydrogen (white). Note the hydrogen bonds between silanol groups on the inner surface (left panel).

  Si-O and Al-O bond lengths, so that the skeleton (aluminum Al, silicon Si, and bridging oxygens O br ) is well reproduced with the CLAYFF potential (in both variants).

Figure 2 :

 2 Figure 2: Left: angle distributions for Si-O int -H int angles (top) and Al-O ext -H ext angles (bottom) for DFT (red line), CLAYFF (dark blue line) and extended-CLAYFF (sky blue line) simulations. Right: density of the internal hydrogen atoms, H int , along the nanotube axis 𝑧 for an empty nanotube (solid lines) and a fully hydrated nanotube (dashed lines). Data from DFT are shown in red, data from CLAYFF in dark blue, and from extended-CLAYFF in sky blue.

Figure 3 :

 3 Figure 3: Fully hydrated imogolite nanotube (𝑛 = 12) with hexagonal packing. Water molecules are represented as blue sticks.

Figure 4 :

 4 Figure 4: Two-dimensional density profiles in the 𝑥𝑦 plane for atoms from the silanols groups (Si, O int , H int ) and water molecules (O w on the left pane, H w on the right pane).

Figure 5 :

 5 Figure 5: Two-dimensional density profiles (the radial dimension versus the projection along the 𝑧 axis) for surface atoms (Si, O int , H int ) and water atoms (O w , H w ). Silicon rings are indicated by dashed lines. Density profiles are normalized, so that for each atom type the highest value is set to 1.0 (full color intensity).

Figure 6 :

 6 Figure 6: Illustration of the model on the 𝑛 = 12 nanotube (only half of the nanotube is shown). The triangular adsorption sites are highlighted by drawing yellow triangles between silicon atoms and red triangles between oxygen atoms. In the right panel, water molecules within the first adsorption layer are drawn and hydrogen bonds are depicted in dotted blue sticks.

Figure 7 :

 7 Figure7: Density profiles on the flattened hydrated nanotube planes for water oxygen, water hydrogen and internal hydrogen. The circular coordinate corresponds to a curvilinear abscissa that draws a circle of radius 𝑅 = 6.5 Å centered on the axis 𝑧. On all the graphs, internal oxygens appear as yellow dots.

Figure 9 :

 9 Figure 9: Illustration of the most common hydrogen bonding patterns for water in an imogolite nanotube. For the numbering of patterns, see Table 2.

Figure 10 :

 10 Figure 10: (a) Autocorrelation of the hydrogen bond existence for three types of donoracceptor couples: water-silanol (black line), silanol-water (blue line) and water-water (red line). (b) Autocorrelation of presence in the first adsorbed layer. The computed data (black line) was fitted by an exponential function (red line), 𝑦 = 1 + 0.42(exp(-𝑥 0.41 ) -1).

Table 1 :

 1 Average number of hydrogen bonds (± the standard deviation) for given donoracceptor pairs along a 1 ns trajectory of empty and hydrated nanotubes. Total number of silanol groups: 72; total number of water molecules: 98.

	Donor Acceptor Empty Nanotube Hydrated Nanotube
	Silanol Silanol	35.4 ± 2.3	2.0 ± 1.3
	Water	Silanol		69.0 ± 2.1
	Silanol	Water		68.7 ± 1.7
	Water	Water		112.4 ± 3.7

Table 2 :

 2 Occurrence of the most recurrent patterns of hydrogen bonding characterized by the number of bonds given to silanol and water and received from silanol and water. Percentages are given with respect to the total number of water molecules.

	Pattern Occurrence Donating Accepting Donating Accepting
			to silanol from silanol to water from water
	1	22.4%	0	0	2	2
	2	14.9%	1	2	1	0
	3	9.3%	2	1	0	1
	4	7.8%	0	0	2	1
	5	7.3%	1	1	1	1
	6	6.3%	1	1	1	0

Table

Table 3 :

 3 Characteristic lifetime for hydrogen bonds calculated from the time length for the autocorrelation function to decay to 50%.

	Donor Acceptor Characteristic time 𝜏
	Water	Water	4.43 ps
	Water	Silanol	57.43 ps
	Silanol	Water	123.5 ps
	the axial diffusion in fully loaded nanotube. In this section, we go further to quantify the
	dynamics of the water inside the imogolite nanotube.

Acknowledgements

We thank Antoine Thill, Sophie Le Caer, and Yuanyuan Liao for fruitful discussions. We acknowledge access to the CINES High Performance Computing platforms provided by a GENCI grant (A0030807069).

Hydrogen bond patterns

The localization of the first layer of adsorbed water on the inner surface of the imogolite nanotube, and the strong hydrogen bonds that are created with the tube's silanol groups, make the first layer of water strongly ordered. We can see on Figure 5 (bottom panel) that this includes rotational ordering of the water molecules, with three marked preferential positions for their hydrogen atoms at room temperature. Two of those, equivalent by symmetry, correspond to the formation of a hydrogen bond, donated by the water molecule to the silanol group (acting as acceptor). The third position corresponds to a hydrogen atom "dangling" toward the inside of the nanotube. Because the water molecules have 2 hydrogen atoms that can occupy these positions, there are 3 possible geometries, summarized on Figure 8: a water molecule can be donating two hydrogen bonds (case A); accepting one and donating one (case B); or accepting two hydrogen bonds from neighboring silanol groups (case C). In this section, we quantified the hydrogen-bond network and dynamics with a set of geometric criteria, that are commonly used in the literature for hydrogen bonds in water.

We consider a hydrogen bond D-H• • • A to be present between donor atom D and acceptor atom A if, and only if, the D• • • A distance is less than 3.5 Å and the ̂ HDA angle is less than 30°. [START_REF] Luzar | Hydrogen-bond kinetics in liquid water[END_REF] We performed hydrogen bond analysis on the hydrated and empty nanotubes, and differentiated the types of hydrogen bonds according to the donor and acceptor groups.

Water hydroxyls and silanols can both be donor and acceptor.

The average number of hydrogen bonds formed during the trajectory is presented in

Supporting Information

Description of the Grand Canonical Monte Carlo adsorption models, adsorption isotherms, additional structural information, and mean-square displacement of adsorbed water.