Stéphane Tesson

Mathieu Salanne

Benjamin Rotenberg

Sami Tazi

Virginie Marry
email: virginie.marry@upmc.fr

A Classical Polarizable Force Field for Clays: Pyrophyllite and Talc

We develop a classical polarizable force field for clays, based on the Polarizable Ion Model (PIM), focusing on two neutral clays: pyrophyllite and talc. The full set of parameters of the force field is determined from density functional theory calculations, using maximally localized Wannier functions with a force-and dipole-optimization procedure. Simulation results for our new polarizable force field (PIM) are compared to the state-of-the-art non-polarizable flexible force field (ClayFF), in order to assess the importance of taking polarization effects into account for the prediction of structural properties. This new force field is validated by comparing the detailed structure of the two minerals against X-ray data. Introducing polarization allows for a good transferability and better agreement with the experimental data, in particular for the subtle orientational deformation of the layers.

Introduction

Clay minerals, 1-3 more precisely phyllosilicates, are minerals that constitute most of the soils and the sedimentary rocks. The wide use of clay minerals in industrial applications (catalytic activities: sorbents, filters, waste deposits), energy and ecological engineering (oil recovery, ground water remediation, geological barrier for radioactive waste and CO 2 ) is partly due to their remarkable properties of retention at the mineral surface. The physico-chemical properties of clay minerals and their interfaces with fluid confined within their pores has been investigated since decades from both the experimental and theoretical points of view.

Classical molecular simulations provide a detailed understanding on the microscopic scale, however their predictions rely on the model used to describe the interactions between the atoms. Therefore, several classical force fields have been developed. [START_REF] Skipper | Computer calculation of water-clay interactions using atomic pair potentials[END_REF][START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF][START_REF] Heinz | Force field for micatype silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF][START_REF] Pitman | Dynamics of confined reactive water in smectite clay-zeolite composites[END_REF] This allowed the use of molecular simulation to describe specific effects, including the prediction of the structure, [START_REF] Delville | Structure and properties of confined liquids: a molecular model of the clay-water interface[END_REF][START_REF] Marry | Microscopic simulations of interlayer structure and dynamics in bihydrated heteroionic montmorillonites[END_REF][START_REF] Wang | Effects of substrate structure and composition on the structure, dynamics, and energetics of water at mineral surfaces: A molecular dynamics modeling study[END_REF][START_REF] Marry | Structure and dynamics of water at a clay surface from molecular dynamics simulation[END_REF] thermodynamics and dynamics of clays, [START_REF] Wang | Effects of substrate structure and composition on the structure, dynamics, and energetics of water at mineral surfaces: A molecular dynamics modeling study[END_REF][START_REF] Marry | Structure and dynamics of water at a clay surface from molecular dynamics simulation[END_REF][START_REF] Tambach | Hysteresis in clay swelling induced by hydrogen bonding: accurate prediction of swelling states[END_REF][START_REF] Botan | How electrostatics influences hydrodynamic boundary conditions: poiseuille and electro-osmostic flows in clay nanopores[END_REF] their hydrophilic/hydrophobic properties [START_REF] Rotenberg | Molecular explanation for why talc surfaces can be both hydrophilic and hydrophobic[END_REF] or the sorption of ions. [START_REF] Greathouse | Molecular dynamics simulation of uranyl(vi) adsorption equilibria onto an external montmorillonite surface[END_REF][START_REF] Greathouse | Water structure and aqueous uranyl(vi) adsorption equilibria onto external surfaces of beidellite, montmorillonite, and pyrophyllite: results from molecular simulations[END_REF][START_REF] Kalinichev | Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: application to cement materials[END_REF] The veracity of results from molecular dynamics simulation strongly depends on the functional form and associated parameters of the force field used to calculate the interactions between atoms. The simulation must correctly describe the balance of interactions between the mineral and the interfacial fluid, reproduce the experimental structure, and capture the effects of ionic sizes and of the polarization of water molecules by ions. The reliability of the force field is even more important for the study of ionic solutions with multivalent ions or at the interface between clays and water molecules, because of multi-body effects. [START_REF] Rotenberg | Molecular explanation for why talc surfaces can be both hydrophilic and hydrophobic[END_REF][START_REF] Salanne | Including many-body effects in models for ionic liquids[END_REF][START_REF] Molina | Ions in solutions: determining their polarizabilities from firstprinciples[END_REF][START_REF] Chang | Recent advances in molecular simulations of ion solvation at liquid interfaces[END_REF] The force fields available in the literature for clays and their interfaces do not take into account the polarizability of molecules, [START_REF] Skipper | Computer calculation of water-clay interactions using atomic pair potentials[END_REF][START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF][START_REF] Pitman | Dynamics of confined reactive water in smectite clay-zeolite composites[END_REF][START_REF] Skipper | Monte Carlo simulation of interlayer molecular structure in swelling clay minerals. I: methodology[END_REF] even though it is known that polarization effects can play a significant role, especially in the presence of an electric field at the interface. [START_REF] Kamath | Comparison of select polarizable and non-polarizable water models in predicting solvation dynamics of water confined between MgO slabss[END_REF] Since the numerical results do not always reproduce quantitatively the experimental data, the relevance of the interpretation of microscopic simulations may be questioned. [START_REF] Marry | Anisotropy of water dynamics in clays: insights from molecular simulations for experimental QENS analysis[END_REF] In the case of pure oxide materials, Madden et al. have developed a series of Polarizable Ion Models (PIM) and Aspherical Ion Models (AIM) from ab initio simulations. [START_REF] Madden | Covalent'effects in 'ionic' systems[END_REF][START_REF] Aguado | Multipoles and interaction potentials in ionic materials from planewave-DFT calculations[END_REF][START_REF] Marrocchelli | The construction of a reliable potential for GeO2 from first principles[END_REF][START_REF] Burbano | A dipole polarizable potential for reduced and doped CeO 2 obtained from first principles[END_REF] In particular, Jahn et al. [START_REF] Jahn | Modeling earth materials from crustal to lower mantle conditions: a transferable set of interaction potentials for the {CMAS} system[END_REF] have determined a full set of AIM parameters for the Ca-Mg-Al-Si-O (CMAS) system, which is the main component of the Earth's crust and mantle. However, this force field cannot readily be extended to aqueous solutions. We have recently shown that polarization effects play an important role in the physico-chemical properties of ions dissolved in water. [START_REF] Tazi | A transferable ab initio based force field for aqueous ions[END_REF][START_REF] Molina | Primitive models of ions in solution from molecular descriptions: a perturbation approach[END_REF][START_REF] Tazi | Description moléculaire des ions aux interfaces argile-eau[END_REF] In the present work, we extend the PIM model to clay minerals. Following the strategy of Madden et al., we derive the parameters of a force field for pyrophyllite and talc.

Contrary to the parametrization of many force fields, the use of Density Functional Theory (DFT) calculations, using Maximally Localized Wannier Functions (MLWFs) [START_REF] Marzari | An introduction to maximally-localized Wannier functions[END_REF] with a forceand dipole-optimization procedure, [START_REF] Salanne | Including many-body effects in models for ionic liquids[END_REF] renders experimental input unnecessary. Clays are layered aluminosilicates, with a large lateral extension compared to their width (∼ 1 nm). We focus here on two specific uncharged clays of the smectite family, namely pyrophyllite [START_REF] Lee | Single crystal X-ray refinement of pyrophyllite-1Tc[END_REF] and talc. [START_REF] Rayner | The crystal structure of talc[END_REF] The manuscript is organized as follows. After introducing the microscopic structures of pyrophyllite and talc, we describe in detail the polarizable force field and its parametrization. Then, the force field is validated by a comprehensive comparison of its prediction of the structure with experimental data.

Structure of pyrophyllite and talc

The Polarizable Force Field

The PIM force field is composed of four different terms: charge-charge, dispersion, overlap repulsion and polarization:

V total = V Charge + V Dispersion + V Repulsion + V P olarization . (1)
The first term corresponds to the electrostatic interaction between two charges,

V Charge = i<j q i q j r ij , (2)
where q i and q j are the charges of each atom and r ij is the distance between them. Here the formal charges are used, they are: O The second term in Eq. 1 corresponds to the dispersion interaction, due to the instantaneous correlations of density fluctuations between the electronic clouds and given by: 35-37

V Dispersion = - i<j f ij 6 (r ij ) C ij 6 (r ij ) 6 + f ij 8 (r ij ) C ij 8 (r ij ) 8 , (3)
where C ij 6 and C ij 8 are the dipole-dipole and dipole-quadrupole dispersion coefficients and r ij is the distance between the two atoms i and j. The Tang-Toennies damping function f ij n is used to correct the short-range interaction as: 38

f ij n (r ij ) = 1 -e -b ij n r ij n k=0 (b ij n r ij ) k k! , (4)
where 1/b ij n is the range of the damping. The third term in Eq. 1 corresponds to the short-range repulsion, described by a simple exponential form:

V Repulsion = i<j A ij e -B ij r ij , (5)
where A ij and B ij are two parameters. Finally, the last term in Eq. 1 correponds to the polarization, composed of 3 contributions: charge-dipole and dipole-dipole interactions, as well as the energy cost for deforming the electronic cloud of the atom:

V P olarization = i<j q i r ij • µ µ µ j r 3 ij g ij 4 (r ij ) - µ µ µ i • r ij q j r 3 ij g ji 4 (r ij ) + µ µ µ i • µ µ µ j r 3 ij - 3 (r ij • µ µ µ i ) (r ij • µ µ µ j ) r 5 ij (6)
+ i |µ µ µ i | 2 2α i ,
where α i is the polarizability of ion i, µ i and µ j are the induced dipole, g ij is the short-range correction to the multipolar expansion by the Tang-Toennies damping function (Eq. 7).

g ij 4 (r ij ) = 1 -c ij e -b ij D r ij 4 k=0 (b ij D r ij ) k k! . (7)
The polarization term includes many-body electrostatic effects since the induced dipoles fluctuate along the simulation depending on the positions of all the ions. They are calculated at each molecular dynamics step by minimizing the polarization energy:

∂V P olarization ∂µ i α = 0. ( 8
)
In the following, we will not go into the details of the two first terms because the parameters of the charge potential and the dispersion potential are already known. [START_REF] Jahn | Modeling earth materials from crustal to lower mantle conditions: a transferable set of interaction potentials for the {CMAS} system[END_REF] We will focus on the two other terms, i.e. repulsion and polarization. These contributions to the overall energy depend on ionic properties such as the ionic radius, for the repulsion, or the condensed-phase polarizabilities. The purpose of the present work is therefore to derive the repulsion and the polarization parameters of the PIM force field for the atomic interactions inside clays. We now describe the procedure used to determine all the parameters from ab initio calculations.

Parametrization DFT calculations

The optimization procedure aims at finding the set of parameters (A ij , B ij , b ij and c ij ) that minimize the error made in the classical calculation of the forces and dipoles with respect to a series of reference DFT calculations. If this error is sufficiently small, the interaction potential can be considered to be of ab initio accuracy. The optimization procedure is composed of three steps.

We first generate a trajectory by classical molecular dynamics simulation with an available force field (see below). This trajectory is set up to obtain some configurations of the clay system. The configurations must be chosen at the thermodynamic equilibrium and separated in order to be independant to each other. Then, ab initio forces and dipoles are computed on a few configurations using DFT-based electronic structure calculations. The ground-state wavefunction {φ • } of each system is obtained by minimization of the Kohn-Sham energy E KS . The force acting on each atom F i , is calculated using:

F i DF T = δE KS [{φ • }] δr i , (9)
where r i is the position of atom i. The Kohn-Sham orbitals are delocalized by nature in the plane-wave basis. This renders the assignment of atomic or molecular properties difficult.

In contrast, the Maximally Localized Wannier Functions (MLWFs) provide a picture of the electronic density around each atom. They are determined by unitary transformations of the Kohn-Sham eigenvectors and are constructed by choosing the phase so that it minimizes their spread. [START_REF] Marzari | An introduction to maximally-localized Wannier functions[END_REF] A complete theory of electric polarization in crystalline dielectrics has been developed in recent years, [START_REF] Souza | Polarization and localization in insulators: generating function approach[END_REF][START_REF] King-Smith | Theory of polarization of crystalline solids[END_REF][START_REF] Vanderbilt | Electric polarization as a bulk quantity and its relation to surface charge[END_REF] which validates the calculation of the dipole moments of single ions or molecules from the center of charge of the subset of MLWFs which are localized in their vicinity. [START_REF] Silvestrelli | Water molecule dipole in the gas and in the liquid phase[END_REF][START_REF] Bernasconi | Cation polarizability from first-principles: Sn2+[END_REF][START_REF] Bernasconi | Ionic to molecular transition in AlCl 3: an examination of the electronic structure[END_REF][START_REF] Souza | Polarization and localization in insulators: generating function approach[END_REF]

Force-and dipole-matching

The final step of the optimization procedure consists in determining the parameters b ij and c ij of the Tang-Toennies function of the polarization term (eq.7) and the parameters A ij and B ij of the repulsion term (Eq. 5) that best reproduce the reference DFT forces and dipoles.

The polarization parameters are determined numerically by minimizing the error function on the dipoles:

χ 2 Dipoles (b ij , c ij ) = 1 N conf 1 N atom conf atom ||µ µ µ classical -µ µ µ DFT || 2 ||µ µ µ DFT || 2 , ( 10
)
where N conf is the number of configurations on which DFT calculations are performed, N atom is the number of atoms per configuration, µ µ µ classical and µ µ µ DF T are the dipoles obtained by classical molecular dynamics using a given set of parameters and by DFT, respectively.

The repulsion parameters are obtained by minimizing:

χ 2 F orces (A ij , B ij , b ij , c ij ) = 1 N conf 1 N atom conf atom ||F classical -F DFT || 2 ||F DFT || 2 , ( 11
)
where F classical and F DF T are the forces obtained by classical molecular dynamics and by DFT calculation, respectively.

The optimization procedure ends when the parameters minimizing χ 2 Dipoles and χ 2 F orces are found. The optimization procedure can be summarized as follows:

1. Generation of a series of representative configurations using Classical MD 2. DFT calculations on each of these configurations (i) Determination of the ground-state wavefunctions, which gives access to the ab initio forces (ii) Wannier localization, from which the ab initio induced dipoles are calculated

Minimization of χ 2

Dipole with respect to the parameters of the polarization term (V P olarization ) and the charge transfer parameter δ within the OH group, and of χ 2 F orce with respect to the repulsion term (V Repulsion ).

Simulations Details

The [START_REF] Lee | Single crystal X-ray refinement of pyrophyllite-1Tc[END_REF][START_REF] Smyth | Crystallographic data for minerals[END_REF] and 9.381 Å for talc. [START_REF] Rayner | The crystal structure of talc[END_REF] Molecular dynamics simulations are performed in the NVT ensemble with the cp2k simulation package 47 at T=300 K. The first configurations are generated using a PIM with an initial choice of parameters. [START_REF] Tazi | Description moléculaire des ions aux interfaces argile-eau[END_REF] During subsequent iterations of the procedure described above, new configurations used for the dipole-and force-matching are generated using a PIM with the current values of these parameters. Periodic boundary conditions are used in the three directions of space. The temperature is controlled via a Nose-Hoover thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular[END_REF][START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF] with a time constant equal to 1 ps. Electrostatic interactions are computed using dipolar Ewald summation, [START_REF] Aguado | Multipoles and interaction potentials in ionic materials from planewave-DFT calculations[END_REF][START_REF] Laino | Notes on "Ewald summation of electrostatic multipole interactions up to quadrupolar level[END_REF] with a tolerance of 10 -7 . For each system, we perform an equilibration of 25 ps followed by a 25 ps production run, using a time step of 0.5 fs in the NVT ensemble.

The parametrization of the force field from ab initio simulations is achieved using N conf = 3 representative configurations of pyrophyllite and talc obtained from the classical molecular dynamics simulations. Density Functional Theory (DFT) calculations are performed on these configurations with the BLYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF] and the PBE 53 functionals for both systems. Goedecker-Teter-Hutter 54-56 pseudopotentials are used with the DZVP and the TZV2P plane-wave basis sets [START_REF] Vandevondele | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases[END_REF] and an energy cutoff of at least 280 Ry. After determining the ground-state wavefunctions, the forces acting on each atom are computed and the dipoles are calculated from the MLWFs. [START_REF] Marzari | An introduction to maximally-localized Wannier functions[END_REF] All ab initio calculations are performed with the cp2k simulation package. [START_REF]CP2[END_REF] The numerical minimization of forces (Eq. 11) and dipoles (Eq. 10) are performed with the Minuit library. [START_REF] James | Minuit -a system for function minimization and analysis of the parameter errors and correlations[END_REF] The final set of parameters depends only slightly on the choice of functional and basis set, as the predicted ab initio dipoles and forces are very similar (see Supplementary Information).

Force field parameters

Figure 5 illustrates the comparison between the dipoles and the forces for one of the pyrophyllite configurations calculated with the classical force field and from the DFT calculations.

Dipoles , is near 17% for pyrophyllite and near 19% for talc. The relative error on the forces, χ 2 F orces , is near 50% for pyrophyllite and near 54% for talc. Note that the overall error is dominated by small dipoles and forces. These results can be considered as a good match, i.e. that the dipoles and the forces on all atoms are satisfyingly reproduced by the classical force field. All the parameters are summarized in Table 2 and3. Table 2: Parameters of the PIM force field for neutral clays. δ is the transferred charge (δ = + 0.8983). The parameters C ij 6 , C ij 8 and b ij n are taken from reference 28. The bond between the hydrogen and the oxygen atoms of the hydroxyl group is rigid.

Systems

Ion pair

A ij B ij C ij 6 C ij 8 b ij n c ij b ij D (ij) (Ha) ( Å-1 ) (Ha Å6 ) (Ha Å8 ) ( Å-1 ) ( Å-1 ) O 2--O

Validation Systems and simulation details

Since the present force field is able to correctly reproduce the ab initio dipoles and forces, we now proceed to its validation against experimental data pertaining to the structure and the density of pyrophyllite and talc. We consider three sizes of pyrophyllite and talc simulation cells to assess the possible influence of finite-size effect. The simulation details are the same as above. Each system is first equilibrated in the NVT ensemble during 50 ps. Then, simulations are performed in the anisotropic NPT ensemble (all box lengths and angles are allowed to evolve independently) under a pressure of 1 bar during 50 ps of equilibration followed by 150 ps of production. The pressure is controled by an extension of the Nose-Hoover barostat developed by Martyna et al [START_REF] Martyna | Nosé-Hoover chains: the canonical ensemble via continuous dynamics[END_REF][START_REF] Martyna | Constant pressure molecular dynamics algorithms[END_REF] with a barostat time constant equal to 9 ps and the same time constant of 1 ps for the thermostat.

Simulations under the same conditions were performed for both our new polarizable force field and the state-of-the-art non-polarizable force field ClayFF. [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF]

Lattice parameters

The unit cell parameters obtained for the equilibrated systems are summarized in Table 5.

Since the effect of system size is very limited, we report here averages over all the systems.

The results obtained with ClayFF or PIM for the parameters of the elementary unit cell for both clays are in good agreement with the experimental data. In particular, both force fields ensure the integrity of the elementary cell and predict the correct lattice parameters and angles. The simulated densities are also compared to the experimental ones in Table 5.

Once again, the overall agreement is good. Such an agreement is already a significant achievement for the present PIM model, which in contrast to ClayFF does not use any experimental input, in particular no information on the experimental structure, during the parametrization process.

Table 5: Unit cell parameters of pyrophyllite and talc clays. a, b and c are the lattice parameters, h is the average interlayer distance, α, β and γ are the angles of the cell and ρ is the density of the system. The relative difference between experimental and simulation results (in percent) is given in parenthesis. The standard errors (SE) are calculated by the error block averaging method. [START_REF] Frenkel | Understanding molecular simulations: from algorithms to applications[END_REF][START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF] The SE, which apply on the last digit of the given values, are estimated to be 3, 8 and 1 for the distances, angles and densities respectively. For example, the 'a' ClayFF parameter of pyrophyllite of 5.216 [3] indicates an SE of 0.003.

Systems Methods a ( Å) b ( Å) c ( Å) h ( Å) α ( • ) β ( • ) γ ( • ) ρ (g.cm -

Sheets structure

Experimentally, the tetrahedral and the octahedral sheets adjust their lateral dimensions in order to accommodate their lateral misfit. As a result, they deviate from an ideal hexagonal symmetry. Figure 6 illutrates the deformation of the hexagonal cavity. The deformation of the structures can be monitored by following well chosen bond lengths and angles. The simulation results for pyrophyllite and talc are summarized in Table S1 and S2, respectively, obtained by both force fields, together with the available experimental data. Before discussing the microstructure in further detail, one can already note the overall good agreement (within a few percents) between simulations with both force fields and with experiments. As mentioned above, it should be noted that the parametrization of the present PIM model

does not rely on experimental information on the microscopic structure. Some discrepancies are observed, which are generally larger with ClayFF (only up to 11-12% in the worst cases). The average thickness of the tetrahedral (∆ Tet. ) and octahedral (∆ Octa. ) of sheets, of the layer (∆ Layer ) and the layer separation (∆ Layer sep. ) obtained for the pyrophyllite and talc equilibrated systems are summarized in Table 6. The results obtained with ClayFF and PIM are both in good agreement with the experimental data. The lateral dimensions of the tetrahedral sheet are usually greater than that of the octahedral one. The distortion of the tetrahedral sheet follows three different mechanisms: [START_REF] Lee | Single crystal X-ray refinement of pyrophyllite-1Tc[END_REF][START_REF] Rayner | The crystal structure of talc[END_REF][START_REF] Martyna | Nosé-Hoover chains: the canonical ensemble via continuous dynamics[END_REF][START_REF] Brigatti | Structures and mineralogy of clay minerals[END_REF] i) tilting of the apical oxygen as evaluated by the angle γ tilt between δ Si-Oa (distance between silicium and apical oxygen atoms) and ∆ Si-Oa (distance between silicium and apical oxygen plans) (Figure 7(a)), ii) rotation of the tetrahedron as evaluated by the angle α rot 3 and4). The relative difference between experimental and simulation results (in percent) is given in parenthesis. The standard errors (SE) are calculated by the error block averaging method. [START_REF] Frenkel | Understanding molecular simulations: from algorithms to applications[END_REF][START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF] The SE, which apply on the last digit of the given values, are estimated to be 7 and 5 for the PIM and ClayFF values respectively. ( -) ( -) ( -) Exp. [START_REF] Rayner | The crystal structure of talc[END_REF] ----

Systems

The tilt angle

The most important distances to consider are the distances between neighbouring apical oxygens (d Oa-Oa ), as illustrated in Figure 6. In fact, a slight extension or shortening of the distance between two apical oxygen atoms will cause the rotation and the tilt of each tetrahedron. Figure 6 illutrates the ideal and the real hexagonal silicate ring arrangements.

In an ideal hexagonal silicate ring, d Oa-Oa = 3.05 Å. For pyrophyllite, [START_REF] Lee | Single crystal X-ray refinement of pyrophyllite-1Tc[END_REF] For both the experimental and the simulation values, the first two distances between 7: Separation between apical oxygen atoms. See figure 6 for the definition of the various oxygen atoms. The standard errors (SE) are calculated by the error block averaging method. [START_REF] Frenkel | Understanding molecular simulations: from algorithms to applications[END_REF][START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF] The SE, which apply on the last digit of the given values, are estimated to be 8 and 7 for the PIM and ClayFF values respectively.

Systems Methods O a (1) to O a (2) via O b (3) ( Å) O a (1) to O a (2) via O b (2) ( Å) O a (1) to O a (2) via O b (1) ( Å)
γ tilt = arccos ∆ Si-Oa δ Si-Oa , (12)
where ∆ Si-Oa is the mean distance between the planes of silicium atoms and that of the apical oxygen atoms, while δ Si-Oa is the mean distance between silicium and apical oxygens. The experimental tilt angle is equal to 6.0 • for pyrophyllite, [START_REF] Lee | Single crystal X-ray refinement of pyrophyllite-1Tc[END_REF] while the PIM and ClayFF force field predict angles of 6.6 • and 4.6 • , respectively (see Table 8). Therefore, including manybody effects via the PIM force field seems to improve the description of this tilt compared to ClayFF. We could not find experimental data on the tilt angle for talc. The tilt angle is smaller in the trioctahedral talc case, because of the cation occupancy of all octahedral sites that hinders the adjustment of the sheet compared to the dioctahedral pyrophyllite case.

The rotation angle

The above-mentioned out-of-plane deformations are also associated with in-plane rearrangements, required to adjust the lateral dimensions of the tetrahedral sheet to that of the octahedral one. This results in a rotation of each tetrahedron by an angle α rot . This rotation decreases the symmetry of the tetrahedral sheet and more specifically from a hexagonal to a di-trigonal silicate ring arrangement (see Figure 6). The mean rotation angle is obtained with the formula:

α rot = |α 1 -120| + |α 2 -120| 4 , (13)
where α rot is the rotation angle, α 8, PIM tends to slightly overestimate the rotation angle by 10-40%, while ClayFF clearly underestimates this distortion of the hexagonal rings, with errors as large as 70-90%. Overall, the treatment of many-body effects within the polarizable PIM force field gives a better account of the flexibility of the mineral structure. This is all the more important that the structural features discussed here are also present in charged clays and that they have a crucial impact on the local structure. [START_REF] Tsipursky | The distribution of octahedral cations in the 2: 1 layers of dioctahedral smectites studied by oblique-texture electron diffraction[END_REF] Moreover, the deformations are known to influence the local properties of the fluid at the surface of the layer, for example the position of the cations or water molecules beside the hexagonal cavities.

Conclusion

We have successfully developed a polarizable force field for clays. The procedure relies only on ab initio DFT calculations; the repulsion and polarization parameters are extracted from Table 8: Average thickness of the silicium-apical oxygen distance ∆ Si-Oa , angles α 1 , α 2 , tilt γ tilt and rotation angle of tetrahedra α rot. (see Figure 7). The relative difference between experimental and simulation results (in percent) is given in parenthesis. The standard errors (SE) are calculated by the error block averaging method. [START_REF] Frenkel | Understanding molecular simulations: from algorithms to applications[END_REF][START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF] The SE, which apply on the last digit of the given values, are estimated to be 8, 5, 6 and 3 for the PIM distances, ClayFF distances, PIM angles and ClayFF angles respectively. The next step in the development of the polarizable force field consists in extending the present approach to the interaction between water molecules, cations (Na + , Cs + , Ca 2+ and Sr 2+ ) and charged clays (montmorillonite and hectorite). Such a force field should provide a more accurate description of these complex charged interfacial systems. In turn, this will improve our understanding of the influence of the deformation on the retention and the dynamics of counter-ions and water molecules at the surface of clays and more generally of the interaction of clay minerals with water, including swelling and wettability properties.

Systems

  Pyrophyllite and talc are two neutral phyllosilicate clays. They are composed of layers piled on top of each other to form particles. Each layer consists of two tetrahedral sheets sandwiching an octahedral sheet (TOT structure). Tetrahedral sheets is composed of basal oxygen (O b ), apical oxygen (O a ) and silicium (Si) atoms, while the octahedral sheet contains apical oxygen atoms, hydroxyl groups (OH) and aluminum or magnesium cations. These sheets are thus made from two cristallographic units: SiO 4 tetrahedra and XO 4 (OH) 2 octahedra (where X is a cation, in our case M = Al or Mg), represented in Figure1.

Figure 1 : 8

 18 Figure 1: Nomenclature of the atoms of the tetrahedral (left) and the octahedral (right) sites.

Figure 2 :

 2 Figure 2: Side and top view of (a) tetrahedral, (b) trioactahedral and (c) dioctahedral sheet.

.

Figure 3 :

 3 Figure 3: Side and top views of the ideal (a) pyrophyllite and (b) talc layers. ∆ Tet. and ∆ Octa. are the average thicknesses of the tetrahedral and octahedral sheets, while ∆ Layer is that of the layer.

Figure 4

 4 illutrates the localization of the Wannier centers in a pyrophyllite configuration. We observe that each oxygen atoms can be attributed four Wannier centers while silicium, aluminium and magnesium atoms do not have any of them in their vicinity. This agrees with previous ab initio MD simulation by Vuilleumier et al.[START_REF] Vuilleumier | Electronic redistribution around oxygen atoms in silicate melts by ab initio molecular dynamics simulation[END_REF]

Figure 4 :

 4 Figure 4: Localization of the Wannier centers in a pyrophyllite configuration. h is the average interlayer distance and ∆ Layer sep. is the average thickness of the layer separation.

Figure 5 :

 5 Figure 5: Dipoles (left) and Forces (right) for each atom for one of the pyrophyllite configurations. The predictions of the classical force field (black lines) for the dipole components (µ x , µ y and µ z ) and the force components (F x , F y and F z ) are compared to the DFT results (red lines).

Figure 6 :

 6 Figure 6: Ideal and real hexagonal silicate ring arrangements

(

  i.e., the deviation from 120 • of each angle in the ring (Figure 7(b) & (c)) and iii) increase in the thickness of the tetrahedral sheet as evaluated by the deviation from 109 • 28' of the O a -Si-O b angle of triads (Figure 7(a)).

Figure 7 :

 7 Figure 7: (a) Definition of the tilt angle of tetrahedra, (b) and (c) bottom-up view of the tetrahedral sheet ; α 1 = 120 + 2 × α rot and α 2 = 120 -2 × α rot

( 1 )

 1 -Si(1)-O b (2)-Si(2)-O a (2) (Figure 7(a)). In turn, the tetrahedra tilt and rotate so that the lateral octahedral and tetrahedral dimensions change. In fact, cations repel each other causing a reconciliation of oxygen atoms common to two octaedra: O a (1)-O a (1) b and O a (2)-O a (2) b (Figure 1). In contrast with the thickening of the tetrahedral sheet, the thickness of the octahedral one decreases. For 2:1 layer silicates, the O a -Si-O b (2) angle is slightly wider than the O a -Si-O b (1) and O a -Si-O b (3) angles. The average of the three O a -Si-O b angles is a measure of the thickness of the tetrahedron: Tetrahedra are slightly elongated along the z direction and compressed in the x and y directions, thereby producing a O a -Si-O b (2) angle larger than the ideal one (109.47 • ). 33 The experimental O a -Si-O b (2) angle is equal to 110.5 • for pyrophyllite 33 and Rayner et al. report a value of 109.2 • for the average O a -Si-O b angle. The average simulated O a -Si-O b (2) angle is equal to 110.5 • and 110.8 • with the PIM and ClayFF force field, respectively, in good agreement with the experimental result. For talc, we obtain 108.9 • and 109.4 • with the PIM and ClayFF force field, respectively, also in good agreement with the experimental value of 109.2 • . As a consequence, the apical oxygen is shifted near the apical oxygens shared by the octahedra edges. This tilt angle γ tilt quantifies the angular deviation of the tetrahedron from the basal plane of apical oxygens (Figure 7(a)). It is given by:

1 and α 2

 2 are respectively the mean angle of O b (2)-O b (3)-O b (1) and O b (1)-O b (2)-O b (3). The PIM (resp. ClayFF) results for the rotation angle are 11.8 • (resp. 1.1 • ) for pyrophyllite and 4.8 • (resp. 0.9 • ) for talc. Compared to the experimental results, reported with the simulation ones in Table

  and dipole-optimization procedure. The force field was built and validated for two neutral clays: pyrophyllite and talc. A full structural description of the atomic clay structure is obtained with the new polarizable ion model (PIM) and compared, together with the stateof-the-art non-polarizable force field for clays (ClayFF), to experimental data. The global structure of the unit cell is well reproduced by both force fields. However, introducing the polarization allows for a better description of the microstructure.

  2-, OH -, Mg 2+ , Al 3+ and Si 4+ .

				Charge transfer within
	the hydroxyl group, of total charge -1, is modelled by partial charges on the corresponding
	atoms: O	(2-δ)-OH	and O	(1-δ)+ H

  The monoclinic talc simulation box contains two clay layers of lateral dimensions 15.87×18.36 Å2 , corresponding to 6 unit cells of formula Si 8 Mg 6 O 20 (OH) 4 per layer. The interlayer spacing is fixed to 9.192 Å for pyrophyllite

orthorhombic pyrophyllite simulation box contains two clay layers of lateral dimensions 20.72 × 26.94 Å2 , corresponding to 12 unit cells of formula Si 8 Al 4 O 20 (OH) 4 per layer.

Table 1 :

 1 χ 2 for the dipoles and the forces for neutral clays

	Systems	χ 2 Dipoles	χ 2 F orces
	Pyrophyllite	0.02981	0.25643
	Talc	0.03609	0.29764

Table 1

 1 shows the error functions, χ 2 Dipoles and χ 2 F orces . The relative error on the dipoles,

	χ

Table 3 :

 3 Electrostatic parameters of the PIM force field for neutral clays.

									2-	28.375 5.486	0.974	5.247	2.718 5.385 4.999
					O 2--Al 3+	60.001 3.326	0.048	0.156	4.168 4.330 2.133
				O 2--Mg 2+	86.277 3.661	0.048	0.156	4.168 4.471 4.996
					O 2--Si 4+	29.402 2.835	0.048	0.156	4.168 5.057 3.165
	Pyrophyllite	O 2--H (1-δ)+	0.000	0.000	0.000	0.000	0.000 2.002 1.861
	and	O (2-δ)-OH	-O 2-	99.000 3.600	0.974	5.247	2.718 3.370 2.439
	Talc	O (2-δ)-OH	-Si 4+	979.802 4.825	0.048	0.156	4.168 3.849 4.269
		O (2-δ)-OH	-Al 3+	20.272 3.146	0.048	0.156	4.168 3.367 1.623
		O	(2-δ)-OH	-Mg 2+	35.623 3.485	0.048	0.156	4.168 3.062 1.455
		O	(2-δ)-OH	-O (2-δ)-OH	79.477 3.300	0.974	5.247	2.718 3.229 4.9823
		O	(2-δ)-OH	-H (1-δ)+	0.000	0.000	0.000	0.000	0.000 2.105 0.124
					Ions	O 2-	OH -
					α i ( Å3 )	0.91	2.39

Table 4

 4

indicates the different box lengths A, B and C and angles α box , β box and γ box .

Table 4 :

 4 Simulated systems for pyrophyllite and talc: A, B and C are the initial box sizes and α box , β box and γ box the initial box angles.

	Systems	Supercell dimensions	A ( Å)	B ( Å)	C ( Å)	α box ( • )	β box ( • )	γ box Number ( • ) of atoms
		5×5×2	25.800 44.830 18.694 90.00 90.00 90.00	2000
	Pyrophyllite	6×6×2	30.960 53.796 18.694 90.00 90.00 90.00	2880
		7×7×2	36.120 62.762 18.694 90.00 90.00 90.00	3920
		5×5×2	26.450 45.865 18.920 90.00 99.00 90.00	2100
	Talc	6×6×2	31.740 55.038 18.920 90.00 99.00 90.00	3024
		7×7×2	37.030 64.211 18.920 90.00 99.00 90.00	4116

Table 6 :

 6 Average thickness of the tetrahedral ∆ Tet. and octahedral ∆ Octa. sheets, layer ∆ Layer and layer separation ∆ Layer sep. (see Figures

  Methods ∆ Tet. ∆ Octa. ∆ Layer ∆ Layer sep.

			( Å)	( Å)	( Å)	( Å)
		PIM	2.15 (0.1) (8.1) 2.25	6.52 ( -)	2.87 (3.9)
	Pyrophyllite	ClayFF	2.04 (5.3) (10.1) ( -) 2.29 6.46	2.86 (3.6)
		Exp. 33	2.15	2.08	-	2.76
		PIM	2.16 ( -)	2.13 ( -)	6.44 ( -)	2.93 ( -)
	Talc	ClayFF	2.00 ( -)	2.32	6.34	2.89

  the mean distance from O a (1) to O a (2) via O b (3) and O a (1) to O a (2) via O b (1), where the definition of these atoms can be found in Figure 6, is 2.78 Å. The third distance, O a (1) to O a (2) via O b (2), is 3.46 Å. For talc, 34 the shortest and largest distances between two apical oxygens are respectively equal to 3.05 Å and 3.31 Å.Table 7 compares the experimental and simulated values. A good agreement is observed for both the ClayFF and PIM force fields.

  shorter than the ideal ones, whereas the third distance is longer. This deformation causes a local deformation of the bridge O a

	apical oxygens are
	ClayFF	2.93	3.29
	Pyrophyllite	PIM	2.81	3.45
		Exp.	2.78	3.46
	ClayFF	3.03	3.11
	Talc	PIM	2.92	3.20
		Exp.	3.05	3.31

Acknowledgement

The authors acknowledge the financial support from the Labex MiChem. We also thank the Centre de Calcul Recherche et Enseignement of Jussieu for access to computational resources.

We are grateful for the computing resources on OCCIGEN (CINES, French National HPC) obtained through the project x2015097321.

Supporting Information Available

Pyrophyllite and talc sheets structures details, comparison of several DFT functional and basis set for the calculation of the reference ab initio dipoles and forces, and the parameters of the ClayFF force field. This material is available free of charge via the Internet at http://pubs.acs.org/.


Stéphane Tesson

Wilfried Louisfrema

Mathieu Salanne

Anne Boutin

Benjamin Rotenberg

Virginie Marry
email: virginie.marry@upmc.fr

Classical Polarizable Force Field to Study Dry Charged Clays and Zeolites

come

Introduction

Clay minerals and zeolites are nanoporous aluminosilicate materials extensively used due to their industrial importance in gas adsorption and separation. Clay minerals [START_REF] Ricq | La classification des minéraux argileux (phyllosilicates)[END_REF][START_REF] Wentworth | A scale of grade and class terms for clastic sediments[END_REF][START_REF]Nickel-Strunz Classification -Primary Groups 10th edition[END_REF] are layered with a large lateral extension compared with their width (∼1 nm). Zeolites 4,5 consist of three-dimensional crystalline frameworks. Both are used in several areas, including health (drugs), domestic products (detergents and water softener), energy and environmental engineering (hydrocarbons cracking for fuel production, retention barrier preventing the release of toxic radioactive species or CO 2 into the biosphere), etc. These applications are made possible thanks to their high specific surface that confer them remarkable physicochemical properties, such as adsorption, retention, cationic exchange, etc. All these properties are modulated by the charge of the mineral framework, compensated by interlayer (clays) or extraframework (zeolites) counterions, and the presence of adsorbed molecules (the location of which depends on the position of compensating cations). Many experimental [START_REF] Radoslovich | Cell dimensions and symmetry of layer-lattice silicates. I. Some structural considerations[END_REF][START_REF] Fitch | Localization of benzene in sodium-y-zeolite by powder neutron diffraction[END_REF][START_REF] Brandani | The effect of water on the adsorption of CO 2 and C 3 H 8 on type x-zeolites[END_REF][START_REF] Pitteloud | The structure of interlayer water in wyoming montmorillonite studied by neutron diffraction with isotopic substitution[END_REF][START_REF] Rinnert | Hydration of a synthetic clay with tetrahedral charges: A multidisciplinary experimental and numerical study[END_REF][START_REF] Ferrage | Investigation of dioctahedral smectite hydration properties by modeling of x-ray diffraction profiles: Influence of layer charge and charge location[END_REF] and computational [START_REF] Angela | Méthodes de simulation moléculaire pour l'étude de la distribution des cations et de l'adsorption de molécules polaires dans les zéolithes[END_REF][START_REF] Kalinichev | Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: Application to cement materials[END_REF][START_REF] Berghout | Density functional theory (DFT) study of the hydratation steps of Na + /Mg 2+ /Ca 2+ /Sr 2+ /Ba 2+ exchanged montmorillonites[END_REF][START_REF] Jeffroy | Molecular simulation of cationic zeolites properties: Thermodynamic and structural properties[END_REF][START_REF] Combariza | Simulating the properties of small pore silicazeolites using interatomic potentials[END_REF][START_REF] Katti | Molecular interactions of kerogen moieties with na-montmorillonite: An experimental and modeling study[END_REF] studies have been performed in order to understand their physicochemical properties. Whereas macroscopic descriptions are abundant (adsorption and exchange isotherm, macroscopic diffusion, etc.), [START_REF] Hirasaki | Wettability: Fundamentals and surface forces[END_REF][START_REF] Fuchs | Adsorption of guest molecules in zeolitic materials: Computational aspects[END_REF][START_REF] Boddenberg | A calorimetric and statistical mechanics study of water adsorption in zeolite NaY[END_REF][START_REF] Jeffroy | Evidence of a framework induced cation redistribution upon water adsorption in cobalt exchanged x faujasite zeolite: A joint experimental and simulation study[END_REF][START_REF] Solc | Wettability of kaolinite (001) surfaces -molecular dynamic study[END_REF][START_REF] Sethia | Adsorption of carbon monoxide, methane and nitrogen on alkaline earth metal ion exchanged zeolite-x: Structure, cation position and adsorption relationship[END_REF] a detailed understanding at the microscopic scale may be difficult to obtain. Nevertheless, combining experimental with molecular simulations can provide insights into the local description of the cationic and molecular adsorption sites, of the environment of counterions, or of the fine structure of the clay and zeolite frameworks.

The reliability of molecular dynamics is based on the correct description of interaction between atoms. In the case of zeolites, most of the available force fields consider the framework as rigid and therefore require prior knowledge of the position of atoms and sites, obtained from experiments. As an example, some of us recently conducted a joint experimental and simulation study of the cation distribution in partially Ni-exchanged Y-faujasite. It was found that the migration of cations upon dehydration is strongly coupled with the framework deformation. This underlines the need for force fields accounting for the framework flexibility. In the case of clays, the force fields available in the literature are able to repro-duce qualitatively the structural properties, [START_REF] Katti | Molecular interactions of kerogen moieties with na-montmorillonite: An experimental and modeling study[END_REF][START_REF] Du | A molecular dynamics simulation study of water structure and adsorption states at talc surfaces[END_REF][START_REF] Marry | Structure and dynamics of water at a clay surface from molecular dynamics simulation[END_REF][START_REF] Rotenberg | Molecular explanation for why talc surfaces can be both hydrophilic and hydrophobic[END_REF][START_REF] Pitman | Dynamics of confined reactive water in smectite clay-zeolite composites[END_REF] the thermodynamics, [START_REF] Kalinichev | Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: Application to cement materials[END_REF][START_REF] Whitley | Free energy, energy, and entropy of swelling in cs-, na-, and sr-montmorillonite clays[END_REF][START_REF] Tambach | Hysteresis in clay swelling induced by hydrogen bonding: Accurate prediction of swelling states[END_REF][START_REF] Mcgrath | Simulating fluid-phase equilibria of water from first principle[END_REF][START_REF] Loganathan | Structure, energetics and dynamics of Cs + and H 2 O in hectorite: Molecular dynamics simulations with an unconstrained substrate surface[END_REF] the dynamics, [START_REF] Marry | Structure and dynamics of water at a clay surface from molecular dynamics simulation[END_REF][START_REF] Pitman | Dynamics of confined reactive water in smectite clay-zeolite composites[END_REF][START_REF] Malikova | Temperature effect in a montmorillonite clay at low hydration -microscopic simulation[END_REF][START_REF] Kumar | swelling, interlayer structure, and hydrogen bonding in organolayered double hydroxides: Insights from molecular dynamics simulation of citrate-intercalated hydrotalcite[END_REF][START_REF] Malikova | Diffusion of water in clays on the microscopic scale: Modeling and experiment[END_REF] the hydrophilic/hydrophobic properties, [START_REF] Rotenberg | Molecular explanation for why talc surfaces can be both hydrophilic and hydrophobic[END_REF] or the sorption of counterions. [START_REF] Kalinichev | Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: Application to cement materials[END_REF][START_REF] Greathouse | Molecular dynamics simulation of uranyl(vi) adsorption equilibria onto an external montmorillonite surface[END_REF][START_REF] Greathouse | Water structure and aqueous Uranyl(vi) adsorption equilibria onto external surfaces of beidellite, montmorillonite, and pyrophyllite: Results from molecular simulations[END_REF] Quantitative agreement with experiments remains however a challenge. Overall, force fields tend to underestimate the interlayer distance [START_REF] Delville | Structure and properties of confined liquids: A molecular model of the clay-water interface[END_REF][START_REF] Marry | Microscopic simulations of interlayer structure and dynamics in bihydrated heteroionic montmorillonites[END_REF][START_REF] Wang | Effects of substrate structure and composition on the structure, dynamics, and energetics of water at mineral surfaces: A molecular dynamics modeling study[END_REF] (in particular for the bivalent couterions), to overestimate the diffusion coefficient of water molecules, [START_REF] Marry | Structure and dynamics of water at a clay surface from molecular dynamics simulation[END_REF][START_REF] Tambach | Hysteresis in clay swelling induced by hydrogen bonding: Accurate prediction of swelling states[END_REF][START_REF] Wang | Effects of substrate structure and composition on the structure, dynamics, and energetics of water at mineral surfaces: A molecular dynamics modeling study[END_REF][START_REF] Bot ¸an | How electrostatics influences hydrodynamic boundary conditions: Poiseuille and electro-osmostic flows in clay nanopores[END_REF] and to overestimate the contact angle between a clay surface and a drop of water. [START_REF] Rotenberg | Molecular explanation for why talc surfaces can be both hydrophilic and hydrophobic[END_REF] The force fields available in the litterature do not take into account the polarizability of molecules, [START_REF] Pitman | Dynamics of confined reactive water in smectite clay-zeolite composites[END_REF][START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF][START_REF] Skipper | Computer calculation of water-clay interactions using atomic pair potentials[END_REF][START_REF] Skipper | Monte Carlo simulation of interlayer molecular structure in swelling clay minerals. I: Methodology[END_REF] which may contribute to the organization of ions and water molecules near the mineral surface where an electric field exists. Such effects were shown to be very important in other fields such as biomolecular simulations. [START_REF] Lopes | Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability: Theory and applications[END_REF][START_REF] Cieplak | Polarization effects in molecular mechanical force fields[END_REF][START_REF] Baker | Polarizable force fields for molecular dynamics simulations of biomolecules[END_REF] We have recently extended the force field based on the Polarizable Ion Model (PIM) to study clay minerals and demonstrated its ability to correctly describe the microstructure of two neutral clays: pyrophyllite and talc. [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] In the present work, we extend the PIM to charged clay minerals. To parametrize the PIM force field, we follow the same strategy described in our previous study, [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] which does not require any experimental input. We focus on different charged clays of the smectite family, namely Na-, Ca-, Sr-and Cs-montmorillonite. Then, we assess the transferability of the PIM force field to other aluminosilicates by considering two different zeolites: MX and MY faujasites (M = Na, Ca and Sr). The manuscript is organized as follows. Sections "Polarizable Ion Model" and "Parametrization of the Force Field" describe the different terms of the polarizable force field and succinctly the parametrization procedure. Then, section "Validation of the Force Field" introduces the microscopic structure of montmorillonites and the force field is validated by comparison between experimental and simulated results. Finally, in Section "Transferability to Zeolites", after a brief introduction of the structure of zeolites, we demonstrate the transferability of the PIM force field to faujasite-type zeolites by studying the cationic distribution.

Polarizable Ion Model

In the PIM model, the potential energy is decomposed into four different terms:

V total = V Charge + V Dispersion + V Repulsion + V Polarization . (1)
The charge term corresponds to the Coulomb interaction (here in atomic units) between two atoms,

V Charge = i<j q i q j r ij , (2)
where q i and q j are the charges of each atom and r ij is the distance between them. In our case, formal charges are used: O .

The dispersion term in Eq. 1 is due to the instantaneous correlations of density fluctuations between the electronic clouds. It is given by: 48-50

V Dispersion = - i<j f ij 6 (r ij ) C ij 6 (r ij ) 6 + f ij 8 (r ij ) C ij 8 (r ij ) 8 , (3)
where C ij 6 and C ij 8 are the dipole-dipole and dipole-quadrupole dispersion coefficients. The Tang-Toennies damping function f ij n is used to correct the short-range interaction as: 51

f ij n (r ij ) = 1 -e -b ij n r ij n k=0 (b ij n r ij ) k k! , (4)
where

1 b ij n
is the range of the damping. The repulsion term in Eq. 1 is modelled using a simple decaying exponential:

V Repulsion = i<j A ij e -B ij r ij , (5)
where A ij and B ij are two parameters. Finally, the polarization term in Eq. 1 is composed of 3 contributions: charge-dipole and dipole-dipole interactions, as well as the energy cost for deforming the electronic cloud of the atom:

V Polarization = i<j q i r ij • µ µ µ j r 3 ij g ij 4 (r ij ) - µ µ µ i • r ij q j r 3 ij g ji 4 (r ij ) + µ µ µ i • µ µ µ j r 3 ij - 3 (r ij • µ µ µ i ) (r ij • µ µ µ j ) r 5 ij (6)
+ i |µ µ µ i | 2 2α i ,
where α i is the polarizability of ion i, µ µ µ i and µ µ µ j are the induced dipoles, g ij is the short-range correction to the multipolar expansion by the Tang-Toennies damping function:

g ij 4 (r ij ) = 1 -c ij e -b ij D r ij 4 k=0 (b ij D r ij ) k k! . ( 7
)
The polarization term includes many-body electrostatic effects since the induced dipoles fluctuate along the simulation depending on the positions of all the ions. They are calculated at each molecular dynamics step by minimizing the polarization energy:

∂V Polarization ∂µ i α = 0. ( 8
)
The purpose of the present work is to derive all the parameters of the PIM repulsion and polarization terms for the atomic interactions between cations and clay layers. The parameters for the interactions between the atoms constituting the sheets and between the cations are taken from preliminary studies. [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF][START_REF] Tazi | Description moléculaire des ions aux interfaces argile-eau[END_REF] We now briefly describe the procedure for obtaining the parameters from ab initio calculations, and refer the reader to our previous work [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] for a complete description.

Parametrization of the Force Field

Optimization Procedure

The optimization procedure aims at finding the set of parameters (A ij , B ij , c ij and b ij D ) that minimize the error made in the classical calculation of the forces and dipoles with respect to a series of reference DFT calculations. To determine A ij and B ij parameters for the repulsion potential (Eq. 5) and c ij and b ij D parameters of the Tang-Toennies function (Eq. 7) for the polarization potential (Eq. 6), we used exactly the same process as for dry pyrophyllite and talc. [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] The optimization procedure can be summarized as follows:

1. Generation of a series of representative configurations using classical MD 2. DFT calculations on each of these configurations (i) Determination of the ground-state wavefunction, which provides the DFT forces (ii) Wannier localization, [START_REF] Silvestrelli | Water molecule dipole in the gas and in the liquid phase[END_REF][START_REF] Bernasconi | Cation polarizability from first-principles: Sn 2+[END_REF][START_REF] Bernasconi | Ionic to molecular transition in AlCl 3 : An examination of the electronic structure[END_REF][START_REF] Souza | Polarization and localization in insulators: Generating function approach[END_REF] from which the DFT induced dipoles are calculated 3. Minimization of the error function on the dipoles χ 2 Dipole with respect to the parameters of the polarization term (V Polarization ) and of χ 2 Force with respect to the repulsion term (V Repulsion ):

χ 2 Dipoles (b ij D , c ij ) = 1 N conf 1 N atom conf atom ||µ µ µ classical -µ µ µ DFT || 2 ||µ µ µ DFT || 2 , ( 9
)
χ 2 Forces (A ij , B ij , b ij D , c ij ) = 1 N conf 1 N atom conf atom ||F classical -F DFT || 2 ||F DFT || 2 , ( 10
)
where N conf is the number of configurations on which DFT calculations are performed, N atom is the number of atoms per configuration, µ µ µ classical are the dipoles obtained by classical molecular dynamics using a given set of parameters. F classical are the forces obtained by classical molecular dynamics.

In this work, we determine the parameters of the repulsion potential (Eq. 5) and of the Tang-Toennies function (Eq. 7) of the polarization potential (Eq. 6) for the interactions between cations and the atoms of the clay layer. All the other interactions, including dispersion terms, are taken from our previous study. [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF]

Simulation Details

The montmorillonite simulation boxes contain two clay layers of lateral dimensions 20.72 The interlayer spacing is fixed to 9.7 Å, 9.5 Å, 9.5 Å and 10.8 Å for Na-montmorillonite, [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites: 2. The Li + , Na + , K + , Rb + , and Cs + exchanged forms[END_REF] Ca-montmorillonite, [START_REF] Cases | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonite. 3. The Mg 2+ , Ca 2+ , Sr 2+ and Ba 2+ exchanged forms[END_REF] Sr-montmorillonite [START_REF] Cases | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonite. 3. The Mg 2+ , Ca 2+ , Sr 2+ and Ba 2+ exchanged forms[END_REF] and Cs-montmorillonite 57 respectively. To generate the initial trajectory (i.e. prior to the force field parameterization), the interactions between cations and basal or hydroxyl oxygen atoms are chosen empirically and taken equal to the interaction between the same cations and the oxygen atom of water from our previous work. [START_REF] Tazi | Description moléculaire des ions aux interfaces argile-eau[END_REF][START_REF] Tazi | A transferable ab initio based force field for aqueous ions[END_REF] The short-range cation-cation interactions are treated by the Lennard-Jones potential, [START_REF] Tazi | Description moléculaire des ions aux interfaces argile-eau[END_REF] but it should be noted that such short-range interactions are not crucial due to the strong electrostatic repulsion between cations. The parameters C ij 6 , C ij 8 et b ij n of the cation-oxygen interactions, of the dipersion potential (Eq. 3), are taken equal to the oxygen-X [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] (where X = Al, Mg or Si). The bond between the hydrogen and the oxygen atoms of the hydroxyl group is rigid. During subsequent iterations of the procedure described above, new configurations used for the dipole-and force-matching are generated using PIM with the current values of these parameters. Molecular dynamics simulations are performed with the version 2.4 of the CP2K simulation package. [START_REF]CP2K developers group[END_REF] Periodic boundary conditions are used in the three directions of space. The temperature T = 300 K is controlled via a Martyna et al. thermostat 61 with a time constant equal to 1 ps. Electrostatic interactions are computed using dipolar Ewald summation, [START_REF] Aguado | Multipoles and interaction potentials in ionic materials from planewave-DFT calculations[END_REF][START_REF] Laino | Notes on "Ewald summation of electrostatic multipole interactions up to quadrupolar level[END_REF] with a tolerance of 10 -7 . For each system, we perform an equilibration of 50 ps followed by a 50 ps production run, using a time step of 0.5 fs.

×
From the classical MD trajectories we sample equilibrium configurations which should be sufficiently separated in order to be independent from each other. The parametrization of the force field is achieved using N conf = 6 representative configurations per montmorillonite (3 tv-Na-montmorillonite and 3 cv-Na-montmorillonite, see below for their descriptions). We have tried to use more configurations when adjusting the parameters, and we observed that it did not influence further the value of the parameters. Density Functional Theory (DFT) calculations are performed on these configurations with the PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] functional for all the systems. Goedecker-Teter-Hutter 65-67 pseudopotentials are used with the DZVP plane-wave basis sets [START_REF] Vandevondele | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases[END_REF] and an energy cutoff of at least 400 Ry. After determining the ground-state wavefunctions, the forces acting on each atom are computed and the dipoles are calculated from the Maximally Localized Wannier Functions 67,69,70 (MLWFs). All these calculations are performed with the CP2K simulation package. [START_REF]CP2K developers group[END_REF] The numerical minimization of forces and dipoles are performed with the Minuit library. [START_REF] James | Minuit -a system for function minimization and analysis of the parameter errors and correlations[END_REF] Force Field Parameters Forces are shown in Table 1. Most of the errors on the forces come from the forces exerted on the aluminum, silicon, magnesium and apical oxygen atoms. All errors observed on the dipoles come from the oxygen atoms of the hydroxyl groups. χ 2 Dipoles and χ 2 Forces obtained for the charged clays are larger than those for the neutral pyrophyllite clay (0.029 and 0.25, respectively [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] ). Most of these relative errors are due to very small values of forces and dipoles. As we will see later, the resulting force field is however able to capture the structural feature of the studied clays.

All the parameters are summarized in Tables 2, 3 and4. The parameters of atoms layer are reported in the supporting information. [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] The relevance of including polarization into the force field to capture many-body effects can be demonstrated by considering the dipoles of different oxygen types in the mineral structure. Specifically, we show in the Supplementary Figure S1 that the dipoles of apical and surface oxygen atoms are very different. Therefore, including the polarizability of these atoms in the force field allows a transferable description of both types of oxygen atoms with the same set of parameters. In a non-polarizable model, these different types of oxygen atoms would typically require different set of parameters.

Moreover, we note that, even in that case, such a description would not allow capturing the change in the charge distribution around surface atoms in the presence of cations, which can be significant in particular with multivalent ions and influences their adsorption at the surface.

Validation of the Force Field

Structure of Montmorillonites

Montmorillonite is a charged clay belonging to the same family as pyrophyllite, namely dioctahedral clays. They are composed of layers piled on top of each other to form particles. Each layer consists of two tetrahedral sheets sandwiching an octahedral sheet (TOT structure).

In neutral pyrophyllite, these sheets are thus made from two crystallographic units: SiO 4 Table 2: Parameters of the PIM force field for the repulsion and polarization terms for the interactions between the cations and the oxygen atoms from the sheet. δ is the transferred charge (δ = + 0.8983).

Damping interaction

Ion pair (ij) in Figures 4b and4c. Depending on the location of the vacancy there are two types of octahedral sheets: 72 i) the two octahedra occupied by Al 3+ cations share an edge with two hydroxyl groups, i.e. occupy so called cis sites, and the remaining trans site is vacant -such an octahedral sheet is thus called trans-vacant (tv-, illustrated in Figures 3c and4b); ii) the two octahedra occupied by Al 3+ cations share an edge with one hydroxyl group and one apical oxygen, i.e. occupy one cis and one trans sites, so that one cis site remains vacantsuch an octahedral sheet is called cis-vacant (cv-, illustrated in Figures 3c and4c). Vantelon et al [START_REF] Vantelon | Iron distribution in the octahedral sheet of dioctahedral smectites. An Fe k-edge x-ray absorption spectroscopy study[END_REF] have shown that substitution occurred predominantly in cis sites for a Wyoming-type montmorillonite. This clay is the one that we have used for our cv--montmorillonite systems, except that we have slightly simplified the system by substituting only the "cis" sites. For this study, we worked with tv-montmorillonite and cv-montmorillonite in order to assess the ability of the polarizable force field to differentiate them. A side view of the montmorillonite systems are given in Figure 5.

A ij B ij C ij 6 C ij 8 b ij n between q i and µ j (Ha) ( Å-1 ) (Ha. Å6 ) (Ha. Å8 ) ( Å-1 ) b ij D c ij ( Å-1 ) (-) O 2--

Systems and Simulation Details

We now proceed to the validation of the force field against experimental data pertaining to the structure of montmorillonite. We consider three sizes of montmorillonite simulation cells for each tv-and cv-clay to assess the possible influence of finite-size effects. lengths and angles are allowed to evolve independently) under a pressure of 1 bar during 150 ps of equilibration followed by 150 ps of production. The pressure is controlled by an extension of the Nose-Hoover barostat developed by Martyna et al [START_REF] Martyna | Nosé-Hoover chains: The canonical ensemble via continuous dynamics[END_REF][START_REF] Martyna | Constant pressure molecular dynamics algorithms[END_REF] with a barostat and a thermostat time constant respectively equal to 2 ps and 1 ps. Simulations under the same conditions were performed for both our new polarizable force field and the state-of-the-art non-polarizable force field ClayFF. [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF] The simulations with the PIM are only around 5 times more computationally expensive than with ClayFF. Further improvements could be obtained by using approaches such as the predictor-corrector algorithm of Kolafa [START_REF] Kolafa | Time-reversible always stable predictor-corrector method for molecular dynamics of polarizable molecules[END_REF] .

Lattice Parameters

The unit cell parameters obtained for the equilibrated systems are summarized in Table 6 and the lattice parameters are exposed on the Figure 6. Since the effect of system size is very limited, we report here averages over all the systems. The results obtained for the elementary unit cell parameters a, b, α, β, and γ with ClayFF or PIM, for all clays, are in good agreement with the experimental data. Both force fields ensure the integrity of the elementary cell.

Experimentally, the distance h ranges between 9.5-9.7 Å for Na-montmorillonite, 57,76-80 9.5-11.65 Å for Ca-montmorillonite, 58,76-79 9.5-10.0 Å for Sr-montmorillonite [START_REF] Cases | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonite. 3. The Mg 2+ , Ca 2+ , Sr 2+ and Ba 2+ exchanged forms[END_REF][START_REF] Norrish | The swelling of montmorillonite[END_REF][START_REF] Ferrage | New insights on the distribution of interlayer water in bi-hydrated smectite from x-ray diffraction profile modeling of 00l reflections[END_REF][START_REF] Ferrage | Investigation of smectite hydration properties by modeling experimental x-ray diffraction patterns: Part I. Montmorillonite hydration properties[END_REF] and 10.7-12.0 Å for Cs-montmorillonite. [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites: 2. The Li + , Na + , K + , Rb + , and Cs + exchanged forms[END_REF][START_REF] Ferrage | Investigation of smectite hydration properties by modeling experimental x-ray diffraction patterns: Part I. Montmorillonite hydration properties[END_REF][START_REF] Rutherford | Effects of exchanged cation on the microporosity of montmorillonite[END_REF] The PIM reproduces the experimental h distance more accurately than the ClayFF for most of cations, i.e. Na + , Ca 2+ and Sr 2+ , while both for force fields predict simular values in good agreement with experiments for Cs + . These different values obtained with PIM and ClayFF may come from the structure of the cavities, as discussed below. Internal Structure of the Layers

The average bond lengths and angles obtained by both force fields for tv-and cv-montmorillonites are summarized in Supplementary Tables S2 andS3, respectively. The average thickness of the tetrahedral (∆ Tet. ) and octahedral (∆ Octa. ) sheets, the layer (∆ Layer ), the layer separation (∆ Layer sep. ), the silicium-apical oxygen distance 47 ∆ Si-Oa , the angles 47 α 1 , α 2 , the tilt γ tilt and the rotation angle 47 of tetrahedra α rot. obtained for the various systems are summarized in Table 7 and shown in Figure 7.

Table 7: Average thickness of the tetrahedral ∆ Tet. and octahedral ∆ Octa. sheets, layer ∆ Layer and layer separation ∆ Layer sep. (see Figure 5), silicium-apical oxygen distance ∆ Si-Oa , angles α 1 , α 2 , tilt γ tilt and rotation angle of tetrahedra α rot. . Standard errors (SE) of 0.07 Å and 1 degres for distances and angles, respectively, are calculated by the block averaging method. [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF][START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF] Force field ∆ Tet. ∆ Octa. ∆ Layer ∆ LayersSep. ∆ Si-Oa γ tilt α 1 α 2 α rot.

( Å) ( Å) ( Å) ( Å) ( Å) ( • ) ( • ) ( • ) ( • ) PIM 2.

Octahedral Sheet

We now examine the differences between tv-and cv-montmorillonite by looking at the structure of the octahedral sheet and investigate the ability of the force fields to capture such OH differs in tv -and cv -montmorillonite. This difference could possibly be due to the orientation of these hydroxyl groups: towards the center of the cav- ity for cv -montmorillonite, in contrast to tv -montmorillonite (Figures 3 and4). In order to minimize the interaction between these two hydroxyl groups, the d 5) Å) whereas it is overestimated in tv -montmorillonite (2.85(4) Å). Thus, the PIM force field better describes the difference between tv -and cv -montmorillonite.

Tetrahedral Sheet

For both ClayFF and PIM, the average bond lengths and angles in Supplementary Tables S2 andS3 are similar to the results obtained on pyrophyllite. [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] Specifically, in contrast to ClayFF the PIM force field is able to capture the deformation of the tetrahedral sheets characterized by the following features: i) the angle O a -Si-O b (2) is slightly wider than the two others O a -Si-O b (1) and O a -Si-O b (3), ii) the tetrahedral sheet adjusts its lateral dimensions, thus causing changes in structure (rotation and tilt of the tetrahedra, which modify the shape of the cavities; this important difference with ClayFF is illustrated in Figure 8) and iii) the average thickness of the tetrahedral layers (∆ Tet. ) is comparable to that in pyrophyllite. Experimentally, on montmorillonites and clays belonging to the same family, Tsipursky et al. [START_REF] Tsipursky | The distribution of octahedral cations in the 2:1 layers of dioctahedral smectites studied by oblique-texture electron diffraction[END_REF] and Drits et al. [START_REF] Drits | Crystal-Chemical factors responsible for the distribution of octahedral cations over trans-and cis-sites in dioctahedral 2:1 layer silicates[END_REF] reported values of ∆ Tet. ) between 2.20 and 2.26 Å, a tilt angle (γ tilt ) between 6.6 and 7.8 • , and a rotation angle of tetrahedra (α rot. ) equal to 11 • . By simulation, as we observed on the neutral clay, [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] the PIM force field predicts values of ∆ Tet. , γ tilt and α rot. in good agreement with experimental ones while the ClayFF force field underestimates these properties. The experimental shape of the cavities is ditrigonal. [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites: 2. The Li + , Na + , K + , Rb + , and Cs + exchanged forms[END_REF][START_REF] Tsipursky | The distribution of octahedral cations in the 2:1 layers of dioctahedral smectites studied by oblique-texture electron diffraction[END_REF][START_REF] Drits | Crystal-Chemical factors responsible for the distribution of octahedral cations over trans-and cis-sites in dioctahedral 2:1 layer silicates[END_REF][START_REF] Beermann | Structure analysis of montmorillonite crystallites by convergent-beam electron diffraction[END_REF][START_REF] Brigatti | Structures and Mineralogy of Clay Minerals[END_REF] It is well reproduced with the PIM force field whereas with the non-polarizable one we obtain hexagonal cavities. These different results are illustrated in Figure 8. This success of the PIM force field is permitted by its better account of the flexibility of the structure. The flexibility also has a crucial impact on the local structure and properties of the interlayer species at the surface of the layer, for example the position of the cations beside the hexagonal cavities.

Cation adsorption

As we previously mentioned, the interlayer distance h is underestimated with the ClayFF force field, whereas it is quite well reproduced with the PIM force field. These differences can now be explained by the shape and dimensions of the cavities and the cations size. Indeed, the shape of the cavities is hexagonal with ClayFF and ditrigonal with PIM. In the hexagonal cavities formed by six basal oxygens, the distance between two opposite basal oxygens is equal to 5.11 (8) Å and the diameter of the inscribed circle (Φ circle in Figure 9) to 4.42(8) Å. The specific shape of ditrigonal cavities results in a smaller space limited by three basal oxygens pointing to the center. The height of this triangular cavity is 3.61 (7) Å and the diameter of the inscribed circle is only 2.41(7) Å. These different sizes of cavities are shown in Figure 9. The cation is located in the center of the inscribed circle (Figure 10). Depending on the shape of the cavity, the interlayer ion may enter more or less deeply (Figure 10). The radii [START_REF] Atkins | Chimie: molécules, matière, métamorphoses[END_REF] of the considered ions are equal to 1.02 Å (Na + ), 1.00 Å (Ca 2+ ), 1.16 Å (Sr 2+ ) and 1.70 Å (Cs + ). The sodium, calcium and strontium cations may enter easily into the hexagonal cavities, whereas in the case of the ditrigonal cavities they are positioned just above. Because cesium is too large, it remains above the cavities in both cases: the results obtained with both force fields are therefore similar for this cation (see Supplementary Figure S2). For all systems, the larger the size of the cation (Na + < Ca 2+ < Sr 2+ < Cs + ) the farther it is from surface. The wider hexagonal cavities predicted by ClayFF can accommodate all cations except for the larger cesium. In contrast, the PIM force field predicts that all ions reside above the narrower di-trigonal cavities, at an increassing distance from the tetrahedral sheet with increasing ionic radius (see Figures S2 andS3). Overall, we can conclude that the better agreement with experiments for the interlayer distance h obtained with PIM compared to ClayFF is due to the better description of the adsorption of the cations inside the cavities, which in turn is permitted by the better account of the framework flexibility. We now investigate the transferability of the PIM to other aluminosilicates, namely zeolites.

Transferability to Zeolites

Unlike two-dimensional phyllosilicates such as clays, zeolites display a three-dimensional framework, consisting only of tetrahedral units, with a diversity of topologies for their assembly. [START_REF]Database of zeolite structures[END_REF] Subsitutions of Si 4+ by Al 3+ ions result as for clays in a negative charge compensated by counterions. Multivalent counterions may deform locally the mineral framework and such deformation impacts the distribution of the cations among the crystallographic sites where they can reside.

Zeolites therefore provide a stringent test on the ability of the PIM force field to correctly fields able to accurately describe the flexibility of the framework, [START_REF] Jeffroy | Molecular simulation of zeolite flexibility[END_REF] since standard force fields typically consider the latter as rigid. [START_REF] Jeffroy | Evidence of a framework induced cation redistribution upon water adsorption in cobalt exchanged x faujasite zeolite: A joint experimental and simulation study[END_REF][START_REF] Louisfrema | Cation redistribution upon dehydration of Na 58 y faujasite zeolite: A joint neutron diffraction and molecular simulation study[END_REF][START_REF] Louisfrema | Cation migration and structural deformations upon dehydration of nickel-exchanged na y-zeolite: A combined neutron diffraction and monte carlo study[END_REF][START_REF] Di Lella | Molecular simulation studies of water physisorption in zeolites[END_REF][START_REF] Abrioux | A molecular simulation study of the distribution of cation in zeolites[END_REF][START_REF] Buttefey | Modélisation et simulation moléculaire des cations extra-charpentes dans les zéolithes de type faujasite[END_REF][START_REF] Jaramillo | New force field for na cations in faujasite-type zeolites[END_REF] This further implies the knowledge of the structure from prior experiments. Specifically, we consider here faujasite (FAU) which is the most widely studied zeolite type. [START_REF]Database of zeolite structures[END_REF] Their crystal structure and the corresponding cationic sites are well known, even though their precise position and occupancy can differ depending on the cations nature and content. [START_REF] Mortier | Compilation of extra frameword sites in zeolites[END_REF]

Simulation Details

The cubic simulation box contains one unit cell (space group Fd3m) of composition M x/n Si 192-x Al x O 384 (see Figure 11), where M n+ refers to sodium (Na + ), calcium (Ca 2+ ), or strontium (Sr 2+ ).

We consider here two different zeolites, Y-faujasite with a Si/Al ratio of 2.3 (x = 58) and X-faujasite with a Si/Al of 1 (x = 96). The distribution of susbstitutions complies with Loewenstein's rule (-Al-O-Al-sequences are prohibited), so that they are alternated in the latter case, but is otherwise random in the former case since it was shown that there is no preferential organization in Y-faujasite. Fig. 1. A stylized drawing of the framework structure of zeolite 3.40 (2) X. Near the center of the each line segment is an oxygen atom.

3.05(2)

The different oxygen atoms are indicated by the numbers 1 to 3.48 (3) 4. Silicon and aluminum alternate at the tetrahedral intersec-2.88 (12) tions, except that Si substitutes for about 4% of the Al. 112.0 (5) Extraframework cation sites are labeled with Roman numerals. 109.9( 5 Sr2+ ions at Sr(2) (site II ) are shown. About 75% of the sodalite 91.5 (3) units have this arrangement. Only three Sr(2) ions surround 69.3(5)/130.2 (5) the remaining 25%. These fractions of 25% and uncertainties 95(2)/162 (2) would likely vanish in this and subsequent captions if this struc-77.8 (7) ture could be extrapolated to Si/Al=1. Ellipsoids of 20% prob-70 (3) ability are shown. Molecular dynamics simulations with the PIM force field are performed using the CP2K package. After a preliminary test, it was found that the parameters describing the interaction of Al 3+ in octahedral units, as obtained in the previous section for clays, were not able to stabilize the zeolite structure in which Al 3+ experiences a tetrahedral environment. For a direct transferability, it would have been necessary to also include such a case (substitutions in the tetrahedral sheet). Rather, here we simply use for aluminum the same parameters as for silicon, except of course for the formal charge. As we now proceed to show, this ansatz, physically motivated by the fact that both ions are isoelectronic, is sufficient to achieve our goals.

In order to properly sample all relevant configurations of counterions at 300 K, which are separated by much larger energy barriers than in the case of clays, we use the simulated annealing method. [START_REF] Aarts | Simulated annealing and boltzmann machines[END_REF][START_REF] Eglese | Simulated annealing: A tool for operational research[END_REF] The cations move significantly away from their initial configurations only at very high temperature. In practice, we perform simulations at 2000 K. However, at such a temperature the framework is not stable, at least as described by the PIM. Therefore, during simulations at 2000 K the framework is kept rigid (NVT simulation, during 50 ps).

Then, the structure is relaxed at a temperature of only 700 K (NPT ensemble at a pressure of 1 bar, during 50 ps). The temperature is then brought back to 300 K where equilibrum properties are sampled, before starting a new annealing cycle (to check the convergence of the cation distribution). The thermostat and barostat time constants are respectively equal to 1 ps and 5 ps.

The extraframework cations may be found in various sites, indicated in Figure 11. In order to assign a given cation to a type of site (I, I', II, II' or III, see below), we use a new method that we introduced recently. [START_REF] Louisfrema | Cation redistribution upon dehydration of Na 58 y faujasite zeolite: A joint neutron diffraction and molecular simulation study[END_REF] Contrary to standard approaches, which requires the knowledge of the precise location of the crystallographic sites from experiments, this method also applies to the case of flexible structures, as it is based on the local environment of the cation (coordination by the different types of oxygen atoms, see Figure 11a).

Framework Structure

The faujasite structure contains three distinct types of cages which are known as supercages, sodalites, and hexagonal prisms (Figure 11). Four distincts oxygen atoms can be found in faujasite. O1 bridging oxygens form the edges of the lateral faces of hexagonal prisms.

O2 and O3 are alternatively found along the edges of the hexagonal faces of these prisms, which connect them to the neighbouring cuboctahedral sodalite cages. More precisely, O3

belong to the square faces of the sodalite cages, while O2 belong to the hexagonal windows between the sodalite and supercages. Lastly, O4 oxygen atoms are located on the shared edges between the square faces and sodalite-supercage hexagonal windows.

As a first test of the performance of the PIM model, we analyze the structure of the mineral framework. More precisely, Tables 8 and9 report the average TO distances (with T=Al or Si) and TOT angles for the various oxygen types, as well as the lattice parameter, for Y and X faujasite, respectively. In each table, the simulation results are compared with experimental data, when available. [START_REF] Fitch | Localization of benzene in sodium-y-zeolite by powder neutron diffraction[END_REF][START_REF] Louisfrema | Cation redistribution upon dehydration of Na 58 y faujasite zeolite: A joint neutron diffraction and molecular simulation study[END_REF][START_REF] Olson | Reinvestigation of the crystal structure of the zeolite hydrated nax[END_REF][100][101] The lattice parameter is in good agreement with experiments, within ∼ 2 % for both Si/Al ratios and all counterions. In addition, the PIM force field is able to capture the lattice contraction with increasing Si/Al ratio (from X to Y). In NaY, all angles and distances are approximately equal (∼ 140

Distribution of Extraframework Cations

Another key feature of zeolites, particularly important for their applications, is the distribution of the extraframework cations among the various sites illustrated in Figure 11a: Site I is at the center of the hexagonal prism; site I' is the image of site I in the sodalite cage; site II' is also inside the sodalite cage, near the center of hexagonal windows; site II is the image of site II' inside the supercage; finally, site III/III' is inside the supercage at/near the center the four-ring placed between two other four-rings.

The results obtained with the PIM force field are reported in Tables 10 and11 for Y and X-faujasites, respectively. These tables also report corresponding experimental data, when available. As can be seen in both tables, the results of the simulations with the PIM model are in excellent agreement with the experimental results for the cation distribution. In NaY almost all sites II (32 sites per unit cell) are occupied by Na + cations, while sites II' and III/III' are generally empty. The remaining cations occupy sites I and I'. In Ca-and Srexchanged Y and X-faujasites, the more confined sites I (and I') are occupied preferentially (16 hexagonal prisms), while the remaining cations reside in sites II.

Finally, we note that in the dry states considered in the present study, the cationic distribution is dominated by interionic repulsion. A more stringent test of the performance of the PIM force field would be the case of hydrated zeolites, which would requires the parametrization of the water-aluminosilicate interaction. We therefore leave this for further investiga-

Conclusion

We have successfully extended the PIM polarizable force field to montmorillonites. We used exactly the same procedure as for neutral clays [START_REF] Tesson | A classical polarizable force field for clays: Pyrophyllite and talc[END_REF] In turn, this will improve our understanding of the interaction of clay minerals with water, including swelling and wettability. Moreover, owing the good performance of the PIM force field for dry Faujasite-type zeolites, we can expect that the parametrization of water-mineral interactions will allow the full study of the simultaneous migration of the cations and the deformation of the framework upon hydration. This approach could also be extended to other aluminosilicates and multivalent ions of industrial and environmental interest (f-block elements, rare Earths, etc.).

Figures 1 and 2

 2 Figures 1 and 2 illustrate the comparison between the dipoles and the forces for one of the montmorillonite configurations calculated with the classical force field and from the DFT calculations. The corresponding error functions, χ 2Dipoles and χ 2 Forces are shown in Table1.

tetrahedra and AlO 4 (Figure 1 :

 41 Figure 1: Forces for each atom for one of the montmorillonite configurations. The predictions of the classical force field (black lines) for the force components (F x , F y and F z ) are compared to the DFT results (red lines).

Figure 2 :

 2 Figure 2: Dipoles for each atom for one of the montmorillonite configurations. The predictions of the classical force field (black lines) for the dipole components (µ x , µ y and µ z ) are compared to the DFT results (red lines).

Figure 3 :

 3 Figure 3: Nomenclature of the atoms. a) tetrahedron, b) octahedron and c) tv-and cv-sites (trans-and cis-vacant respectively).

Figure 4 :

 4 Figure 4: Top and side views of a) tetrahedral, b) tv-octahedral and c) cv-octahedral sheet.

Figure 5 :

 5 Figure 5: Side view of the ideal tv-montmorillonite layers. h is the average interlayer distance, ∆ Layer sep. is the average thickness of the layer separation, ∆ Tet. and ∆ Octa. are the average thicknesses of the tetrahedral and octahedral sheets, while ∆ Layer is that of the layer.

Figure 6 :

 6 Figure 6: Lattice parameters.

  differences. More precisely, we analyze the distance d O OH -O b OH between the oxygens of the hydroxyl groups of the same edge of an octahedron (Figure 3). Experimentally, Wilson et al. have shown that d O OH -O b

Figure 7 :

 7 Figure 7: (a) Definition of the tilt angle of tetrahedra, (b) and (c) bottom-up view of the tetrahedral sheet ; α 1 = 120 + 2 × α rot and α 2 = 120 -2 × α rot .

Figure 8 :

 8 Figure 8: Ring arrangements of montmorillonites obtain with the ClayFF et the PIM force field.

Figure 9 :

 9 Figure 9: Form and size of the hexagonal and ditrigonal cavity. The blue dotted lines represent the periphery of the cavity.

Figure 10 :

 10 Figure 10: Top and side view of the position of the Na cation above the cavity.

Fig. 2 .

 2 A stereoview of a sodalite cavity with attached D6R of ) 88.5(3)/91.5(3) dehydrated Sr 46 -X. One Sr2+ ion at Sr(1) (site I ) and four 102.6(3)

Figure 11 :

 11 Figure 11: (a) Framework structure of faujasite zeolites. 96 There are 4 nonequivalent oxygen atom types, labelled from 1 to 4 (see text). Extraframework cationic sites are labelled I, I', II, II' and III. (b) Snapshot of faujasite structure used for molecular simulations.

1 ± 1 24

 11 • and ∼ 1.7 Å, respectively) and in good agreement with experiments. The bond lengths are slightly shorter, in CaY and SrY compared to NaY.In addition, the tetrahedra are now distorted, with different TOT angles depending on the considered oxygen type. In particular, the TO(1)T, TO(2)T and TO(3)T angles decrease, while the TO(4)T angle slightly increases from NaY to Ca/SrY, in line with the experimental results. This local deformation of the structure may be at the origin of the larger standard deviation on the bond lengths (compared to NaY), as well as of the 1-2 % difference in the lattice parameter with respect to experiments. Consistently with the evolution of the lattice parameter with Si/Al ratio, one observes a slight shortening of the TO bond lengths from CaX to CaY and from SrX to SrY. Overall, the PIM model which was not parametrized on zeolites is able to correctly describe the structure of the faujasite framework, including its evolution with Si/Al ratio (at least in the considered range, i.e. 1-2.3) and its deformation induced by divalent cations. This is in itself a significant achievement of the present PIM. 0.7 137.5 ± 0.7 136.0 ± 0.8 137.5 ± 0.6 1.69 ± 0.05 1.73 ± 0.1 1.73 ± 0.1 1.73 ± 0.

  in order to obtain the force field parameters. The PIM force field was built and validated for two charged clays: tv-and cv-montmorillonites with four different cations (Na + , Ca 2+ , Sr 2+ and Cs + ). A structural description of the atomic clay structure is obtained with the new polarizable ion model (PIM) and compared, together with the state-of-the-art non-polarizable force field for clays (ClayFF), to experimental data. The global structure of the unit cell is well reproduced by both force fields. However, introducing the polarization allows for a better description of the structural flexibility and the physisorption of cations at the tetrahedral sheet surface.Moreover, we have demonstrated the transferability of the PIM model for other aluminosilicates, at the price of a slight adjustement of the force field, namely using the Si-O short range interactions for tetrahedral Al, instead of the parameters determined for octahedral Al in montmorillonite clay. To this end, we considered the most widely studied zeolites, of faujasite-type. Specifically, we investigated NaY and completely exchanged X and Y faujasites with Ca 2+ and Sr 2+ cations. A good agreement with experiments was obtained for the mineral framework structure, including the variations with Si/Al ratio and change in extraframework cation, which is already an achievement in itself. Furthermore, the predicted distribution of extraframework cations is also in excellent agreement with experiments. The next step in the development of the polarizable force field consists in extending the present approach to hydrated Na-, Cs-, Ca-and Sr-montmorillonites. The PIM force field should then provide a more accurate description of liquid-solid interfaces in such systems.

  17.96 Å2 , corresponding to 8 unit cells of formula X 0.75/n Si 8 Al 3.25 Mg 0.75 O 20 (OH) 4 per layer (where X = Na + , Ca 2+ , Sr 2+ or Cs + and n is the formal charge of counterions).

Table 1 :

 1 χ 2 for the dipoles and the forces for charged clays

	Systems	χ 2 Dipoles	χ 2 Forces
	Na-montmorillonite	0.229	0.787
	Ca-montmorillonite	0.073	0.815
	Sr-montmorillonite	0.061	0.790
	Cs-montmorillonite	0.076	0.887

Table 3 :

 3 Polarization damping parameters of the PIM force field. δ is the transferred charge (δ = + 0.8983).

		Na +	30.002	3.300	0.048	0.156	4.168	4.330	2.133
	O 2--Ca 2+	20.004	2.543	0.048	0.156	4.168	5.999	0.028
	O 2--Sr 2+	20.064	2.520	0.048	0.156	4.168	5.769	2.178
	O 2--Cs +	20.001	2.748	0.048	0.156	4.168	5.943	0.219
	O	(2-δ)-OH	-Na + 4019.112 5.998	0.048	0.156	4.168	5.678	3.468
	O	(2-δ)-OH	-Ca 2+ 4998.674 5.654	0.048	0.156	4.168	5.999	0.001
	O	(2-δ)-OH	-Sr 2+ 4994.966 5.419	0.048	0.156	4.168	5.992	0.001
	O	(2-δ)-OH	-Cs +	4000.000 5.100	0.048	0.156	4.168	5.686	0.006

Table 4 :

 4 Atomic polarizabilities.

	Ions	O 2-	O -OH	Na +	Ca 2+	Sr 2+	Cs +
	α i ( Å3 )	0.91	2.39	0.18	0.44	0.81	2.02

hedral sheet contains apical oxygen and aluminium (Al) atoms, as well as hydroxyl groups (OH). In montmorillonite, some octahedral Al are replaced by magnesium (Mg), resulting in a negative charge of the sheet compensated by counterions located in the interlayer space (here Na + , Ca 2+ , Sr 2+ or Cs + ).

Table 5

 5

indicates the different box lengths A, B and C. The simulation details are the same as above. The simulations are first performed in the NVT ensemble during 100 ps in order to equilibrate each system. Then, simulations are performed in the anisotropic NPT ensemble (all box

Table 5 :

 5 Simulated systems for montmorillonite: A, B and C are the initial box sizes.

	Systems	Supercell	A	B	C	Number
	tv-and cv-	dimensions ( Å)	( Å)	( Å) of atoms
		8×4×2		35.88		2608
	Na-montmorillonite	8×5×2	41.28 44.85 19.40	3260
		8×6×2		53.82		3912
		8×4×2		35.88		2584
	Ca-montmorillonite	8×5×2	41.28 44.85 19.00	3230
		8×6×2		53.82		3876
		8×4×2		35.88		2584
	Sr-montmorillonite	8×5×2	41.28 44.85 19.00	3230
		8×6×2		53.82		3876
		8×4×2		35.88		2608
	Cs-montmorillonite	8×5×2	41.28 44.85 21.60	3260
		8×6×2		53.82		3912

Table 6 :

 6 Unit cell parameters of montmorillonite (averaged over tv-and cv-montmorillonite and over the various supercells): lattice parameters a, b, c and cell angles α, β, γ. Standard errors (SE) of 0.04 Å and 0.8 degres, respectively, are calculated by the block averaging method.[START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF][START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF]

		Na +	5.19	9.00	9.77	9.61	89.9	98.9	90.1
	PIM	Ca 2+	5.21	8.97 10.08 10.00	88.8	97.2	90.2
		Sr 2+	5.22	9.01 10.15 10.04	90.0	98.3	90.1
		Cs +	5.22	8.96 10.94 10.69	89.9	98.5	90.0
		Na +	5.20	9.00	9.54	9.41	88.5	97.2	90.0
	ClayFF	Ca 2+	5.21	9.04	9.59	9.41	90.8	98.6	89.8
		Sr 2+	5.20	9.04	9.82	9.59	90.8	99.0	90.1
		Cs +	5.19	9.00 10.98 10.75	90.0	99.5	90.0
	Exp. 72	Na +	5.18	8.98			90.0	99.5	90.0
	Exp. 57,76-80	Na +				9.5-9.7

Force field

Counterion a ( Å) b ( Å) c ( Å) h ( Å) α ( • ) β ( • ) γ ( • )

  O OH -O b OH distance is larger in cv -montmorillonite (2.85-2.88 Å) than in the tv -montmorillonite (2.40-2.50 Å). By simulation, the distances d O OH -O b OH are well reproduced for both clay types with the PIM force field, with d O OH -O b OH = 2.89(5) and 2.42(4) Å for cv -and tv -montmorillonite, respectively. In the case of ClayFF, the distance d O OH -O b OH is reproduced in cv -montmorillonite (2.92(

Table 8 :

 8 Interatomic angles and distances, and lattice parameter for dehydrated Y-faujasite zeolites. Simulation results with the PIM force field are compared with experimental data, when available.

	24.85 ± 0.03
	1.65
	1.67

Table 9 :

 9 Interatomic angles and distances, and lattice parameter for dehydrated X-faujasite zeolites. Simulation results with the PIM force field are compared with experimental data, when available.

Table 10 :

 10 tion. A recent study of cation migration upon (de)hydration in partially Ni-exchangedNa-faujasite 90 demonstrated the importance of accounting for the framework deformation to properly predict the distribution of cations among the various sites and its evolution with water content. Even though the purpose of the present work was to demonstrate the trans-Distribution of cations between the various sites (number of cations per unit cell) for dehydrated Y-faujasite.Simulation results with the PIM force field are compared with experimental data, when available.

	ferability of the PIM force field, the results are also promising for future studies dedicated
	to zeolites.

Table 11 :

 11 Distribution of cations between the various sites (number of cations per unit cell) for dehydrated Y-faujasite. Simulation results with the PIM force field are compared with experimental data, when available.

Acknowledgement S.T. acknowledges financial support from Labex MiChem, W.L. from Région Ile-de-France via the DIM OXYMORE and V.M. from NEEDS-MIPOR via project TRANSREAC. We also thank the Direction des Systèmes d'Information of Jussieu for access to computational resources. We are grateful for the computing resources on OCCIGEN (CINES, French National HPC) obtained through the projects x2015097321 and x2016087684.

Supporting Information Available

Parameters of the PIM force field obtained on pyrophyllite and talc, Trans-and Cis-vacant Na-, Ca-, Sr-and Cs-montmorillonite sheets structures details, Cations adsorption on the tetrahedral sheet surface, and the parameters of the ClayFF force field. This material is available free of charge via the Internet at http://pubs.acs.org/.