Classical Polarizable Force Field to Study Hydrated Hectorite: Optimization on DFT Calculations and Validation against XRD Data

: Following our previous works on dioctahedral clays, we extend the classical Polarizable Ion Model (PIM) to trioctahedral clays, by considering dry Na, Cs-, Ca-and Sr-hectorites as well as hydrated Na-hectorite. The parameters of the force ﬁeld are determined by optimizing the atomic forces and dipoles on density functional theory calculations. The simulation results are validated by comparison with experimental X-ray diffraction (XRD) data. The XRD patterns calculated from classical molecular dynamics simulations performed with the PIM force ﬁeld are in very good agreement with experimental results. In the bihydrated state, the less structured electronic density proﬁle obtained with PIM compared to the one from the state-of-the-art non-polarizable force ﬁeld clayFF explains the slightly better agreement between the PIM results and experiments.

However, they do not account for the polarizability of the atoms and molecules, which can play an important role at the interface between the fluid and the charged layers and especially if the solution contains multivalent and/or ions with large radii [28][29][30].We expect that taking the polarizability into account could provide new quantitative insights into the behavior of clays in the presence of water.In particular, the study of synthetic fluorohectorites with very well defined swelling behavior compared to natural clays (no interstratification) showed low water uptakes and diffusion coefficients which could not be reproduced by molecular dynamics simulations with existing force fields [9,31,32].This suggests that the behavior of the fluid in clay nanopores could be rather different from the one obtained from molecular simulations and the role of the polarizability on the dynamics in clays was already pointed out in the literature [28].We have recently extended the force field based on the Polarizable Ion Model (PIM) [33,34] to two neutral clays, talc and pyrophyllite, and a charged equivalent of the latter, a montmorillonite compensated by cations of various charge and size (Na + , Cs + , Ca 2+ and Sr 2+ ) [35,36].The simulations performed with the PIM model were shown to give results equivalent or in better agreement with X-ray diffraction data than the ones performed with the popular and widely used force-field clayFF.In particular, the deformation of the clay layer and the layer-to-layer distances were better retrieved.Moreover we have shown that the PIM model could be transferred to other types of aluminosilicates such as zeolites [36].
In the present article, we focus on another type of clay mineral, an hydroxylated hectorite (OH-Ht).The studied hectorite is the charged equivalent of talc, with Na + , Cs + , Ca 2+ and Sr 2+ as counter-ions.The first results obtained in the presence of water are also presented.Talc and hectorite have a TOT structure: they are composed of two tetrahedral sheets of SiO 4 sandwiching an octahedral sheet of XO 4 (OH) 2 where X is a cation, here a magnesium Mg 2+ .The elementary cell of talc is Si 8 Mg 6 O 20 (OH) 4 .In hectorite, a fraction of Mg 2+ is substituted by cations of lower charge, here Li + .Thus, the elementary cell of hectorite is X x/n Si 8 Mg 6−x Li x O 20 (OH) 4 [31].X is a cation of valence n which compensates the negative charge of the clay sheet.
In talc and hectorite (trioctahedral clays), all the octahedral sites present in the octahedral layers are occupied by a magnesium or a lithium whereas in pyrophyllite and montmorillonite (dioctahedral clays), only two thirds of the octahedral sites are occupied.As a consequence, the orientation of hydroxyl groups of the octahedral layer are different in these two types of clays.In the case of talc and hectorite, hydroxyl groups have their hydrogen atom pointing away from the surface of the clay layer as seen in Figure 1, which gives the possibility for water molecules to adsorb at low humidity by forming hydrogen bonds with the hydroxyl [6]: these hydrogen bonds are responsible for the hydrophilic property of talc at low relative humidity [37] which is not seen in the case of pyrophyllite or fuorotalc, which remain hydrophobic at all humidities [6,37].In charged clays however, the ability of the cations to hydrate has a crucial influence on the interaction with water and explains why hectorites are hydrophilic and have similar swelling behaviours to other smectites such as montmorillonite [31,38].In the present article, we determine the PIM parameters for an hectorite with a charge close to x = 0.7.In the dry state, several counter-ions are considered: Na + , Cs + , Ca 2+ and Sr 2+ .In the hydrated state, we focus on the hectorite containing one or two layers of water (mono-and bihydrated states) with Na + as counterions.The PIM force field for the hydrated systems is then validated by comparison with X-ray diffraction (XRD) patterns.Indeed, combinations of molecular simulations and diffraction methods have been shown to provide constraints on the computed interlayer model and associated force field used to model clay-water systems [39].Since both methods probe the periodic structures on the same length scale, their combination has allowed a deeper understanding of interlayer organization of water in smectites with different compositions or hydration states [7,40,41] as well as in the presence of organic molecules in the interlayer [42,43].There exists relatively few computational studies of hydrated hectorites [10,[44][45][46][47][48][49].Our results were also compared with the simulation results given by the state-of-the-art force field ClayFF [27] and previous relevant studies.

Characteristics of the Polarizable Ion Model
In the PIM model, the potential energy is decomposed into four different terms: ( The charge term corresponds to the Coulomb interaction (here in atomic units) between two atoms, where q i and q j are the charges of the atoms i and j respectively, and r ij is the distance between them.In PIM model, the formal charges are used: O with δ = +0.8983.The dispersion term in Equation ( 1) is due to the instantaneous correlations of density fluctuations between the electronic clouds.It is given by [50][51][52]: 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 [53]: where 1/b ij n is the range of the damping.The repulsion term in Equation ( 1) is modelled as: where A ij and B ij are two parameters.Finally, the polarization term in Equation (1) contains three contributions: charge-dipole and dipole-dipole interactions, as well as the energy cost for deforming the electronic cloud of the atom: 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: 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 dynamic step by minimizing the polarization energy: The purpose of the present work is to derive the parameters of the PIM repulsion and polarization terms for the atomic interactions between cations, water and clay layers.Some of these parameters have already been determined in our previous studies: the interactions between the counter-ions and the atoms constituting the sheets of montmorillonites [35,36] and the interactions between Na + and water [54] are already known.The interactions parameters between water and clay have been determined in hydrated montmorillonites and constitute the subject of an article in preparation [55].Moreover, the parameters for the dispersion term C ij 6 , C ij 8 and b ij n were taken from elsewhere [34,35].Repulsive and dispersion interaction between cations were neglected due to the predominant strong electrostatic repulsion.In this study we focus on the repulsive and polarization interactions between the lithium present in the octahedral layers of hectorites and all the other atoms of the system.In the following, we briefly describe the procedure for obtaining the parameters from ab initio calculations, and refer the reader to our previous work [35] for a complete description.

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 Density Functional Theory (DFT) calculations.To determine A ij and B ij parameters for the repulsion potential (Equation ( 5)) and c ij and b ij D parameters of the Tang-Toennies function (Equation ( 7)) for the polarization potential (Equation ( 6)), we used exactly the same process that was used previously [35,36].The optimization procedure can be summarized as follows: 1. Generation of a series of representative configurations using classical molecular dynamics (MD) 2. DFT calculations on each of these configurations (i) Determination of the ground-state wavefunction, which provides the DFT forces (ii) Wannier localization [56][57][58][59], from which the DFT induced dipoles are calculated

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 ): 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.

Dry Hectorites
Four hectorite boxes were used for the optimization of PIM parameters.For systems with monovalent counter-ions, they contained two clay layers of lateral dimensions 26.30Å × 18.20 Å each, corresponding to ten unit cells of formula X 0.7 Si 8 Mg 5.3 Li 0.7 O 20 (OH) 4 with X = Na + or Cs + .For systems with divalent counter-ions, they contained two clay layers of lateral dimensions 21.04 Å × 18.20 Å each, corresponding to eight unit cells of formula X 0.75/2 Si 8 Mg 5.25 Li 0.75 O 20 (OH) 4 with X = Ca 2+ or Sr 2+ .The lithium cations were placed randomly in the octahedral sheet, with an exclusion rule preventing two substitutions from being adjacent.During the optimization, the layer-to-layer distances were fixed to 10 Å, typical values of dry layer-to-layer distances of smectites being situated between 9.5 and 10.5 Å [40,[60][61][62][63] hectorite being no exception with a distance of about 9.7 Å for sodium cations [31,64].To generate the initial trajectories from which the configurations used for the optimization are taken, the parameters for lithium interactions with clay oxygen atoms were chosen to be the same as between a hydrated Li + and water oxygen atoms [54].The starting parameters between polarizable cations Ca 2+ , Sr 2+ , Cs + , Na + and the lithium charge, c Li−Cation and b Li−Cation D , were arbitrarily set to zero.The polarizability of Li was neglected [54].All the other parameters were determined previously and can be found in Supplementary Materials.
Molecular dynamics (MD) simulations were performed with version 2.4 of CP2K simulation package [65].Periodic boundary conditions were used in the three directions of space.The temperature T = 300 K was controlled via a Martyna et al. thermostat [66] with a time constant equal to 1 ps.Electrostatic interactions were computed using dipolar Ewald summation [33,67], with a tolerance of 10 −7 .For each system we performed an equilibration of 50 ps using a time step of 0.5 fs.The last configuration was taken for the optimization.Thus, the optimization of lithium parameters was done on four configurations at the same time, each configuration containing a different counter-ion.We have shown previously that three configurations (with a number of atoms comparable to the present case) are sufficient to ensure convergence and that increasing the number of configurations does not change the parameters.In the present case, the fact that the counter-ions are different in the four configurations is not critical since the interactions which are crucial for describing Li + behavior are the repulsion (Equation ( 5)) and damping terms (Equation ( 7)) between Li + which is located in the middle of the sheet and its surrounding clay oxide anions.
Density functional theory (DFT) calculations were performed on these configurations with the PBE functional for all of the systems [68].Goedecker-Teter-Hutter pseudopotentials [69][70][71] were used with DZVP plane-wave basis sets and an energy cutoff of at least 400 Ry.After determining the ground state wave functions, the forces acting on each atom were computed, and the dipoles were calculated from the Maximally Localized Wannier Functions (MLWFs) [71][72][73].The numerical minimization of forces and dipoles were performed with the Minuit library [74].

Hydrated Na-Hectorite
For hydrated hectorites, we used the rigid Dang-Chang polarizable model for water [75].The Dang-Chang model is a four-site model, the fourth site MW carrying the negative charge and the dipole of the molecule.The characteristics of the model are given in Table 1.The PIM parameters for the hydrated cations have been determined previously and have shown to give excellent agreement with the experiment from the structural, thermodynamical and dynamical points of view [54].This model has also been used to parametrize the PIM model in montmorillonite [55].The SHAKE algorithm was used to enforce the rigidity of the water molecules [76].
The optimization of lithium parameters was done in the bihydrated state of the Na-hectorite, i.e., the interlayer spaces contained two layers of water.The charge and the lateral dimensions of the clay layers were the same as in the dry systems.The layer-to-layer distances were fixed to the experimental value of 15.51 Å corresponding to an hectorite of similar charge at 80% of humidity [31].The number of water molecules per unit cell was fixed to the experimental value of 8.9, which corresponds to 178 water molecules in the box.To generate the initial trajectories from which the configurations used for the optimization were taken, the parameters for lithium interactions with water oxygen atoms were chosen to be the same as between a hydrated Li + in solution and water oxygen atoms [54].

Force Field Parameters
Figure 2 illustrates the comparison between the dipoles and the forces for one of the hydrated hectorite configurations calculated with the optimized classical force field and from the DFT calculations.
The corresponding error functions for the dehydrated and the bihydrated hectorite, χ 2 Dipoles and χ 2  Forces are shown in Table 2. Most of the errors on the forces come from the forces exerted on the aluminium, 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 with the dry hectorites are larger than those with the neutral clays [35] but similar to the errors found for dry montmorillonites from our previous work [36].The errors on hydrated charged clays are larger, especially for the forces.χ 2 Dipoles and χ 2 Forces coming from the fit of the interactions between water and clay in hydrated montmorillonite were found to be 0.073 and 1.253 respectively [55].These high relative errors are due to very small values of forces and dipoles.For hectorite, all the parameters were taken to be the same as in montmorillonite, except for Li.Indeed, we expect PIM to be transferable from one clay to another, as long as the environment of the atoms remains the same, which is not exactly the case since Mg 2+ represent only a part of the metallic cations in montmorillonite, whereas they are the majority in hectorite.As a consequence, keeping the parameters obtained on montmorillonite for hectorite leads to higher errors for hectorite.As shown below, the resulting force field is however able to capture the structural characteristics of the studied clays.All the parameters obtained for Li are summarized in Table 3.   [34].The repulsion and dispersion interactions between Li and the other cations are negligible because the electrostatic repulsion predominates.As lithium polarizability is negligible, the polarization potential is restricted to the interactions between the charge of Li and the dipoles of polarizable ions.

Validation of the Force Field
In order to proceed to the validation of the force field, molecular dynamics simulations were performed with the parameters determined above.The validation was done by comparison with existing XRD data.For the hydrated systems, MD-generated XRD profiles were compared to experimental XRD patterns [31] obtained on a hydroxylated hectorite of similar charge.

Simulations Details
In order to obtain the characteristics of the unit cells in dry hectorites, simulations were performed with CP2K package [65] in the anisotropic NPT ensemble, after a phase of equilibration in the NVT ensemble.The pressure was fixed to 1 bar thanks to an extension of the Nosé-Hoover barostat developed by Martyna et al. [66] with barostat and thermostat time constants equal to 2 and 1 ps respectively.The other simulation details are the same as the ones used for the generation of configurations described above.
For hydrated hectorites, 1 ns production runs were performed in the NVT ensemble by fixing the layer-to-layer distances and the number of water molecules to the experimental values obtained by Dazas et al. [31] and corresponding to the conditions of relative humidity detailed in the experimental section below.The characteristics of the simulation boxes are given in Table 4.For comparison, simulations were done on the same systems with the state-of-the-art nonpolarizable force field clayFF [27] and the TIP4P/2005 water model [77] which is one of the best models to describe the structural and dynamical behavior of water on a large range of pressure and temperature.The simulations with clayFF were performed with the LAMMPS package [78].Periodic boundary conditions were used in the three directions of space.The temperature T = 300 K was controlled via a Nosé-Hoover style thermostat [79] with a time constant equal to 1 ps.Electrostatic interactions were computed using particle-particle particle-mesh solver [80], with an accuracy of 10 −5 .The SHAKE algorithm was used to enforce the rigidity of the water molecules [76].A realistic configuration was first generated with a home made Monte Carlo code in which clay layers are rigid.Then 5 ns runs were used for the analysis of the trajectories after 50 ps of equilibration.

XRD Materials and Methods
The hectorite sample with structural formula Na 0.72 Si 8 Mg 5.28 Li 0.72 O 20 (OH) 4 was synthesized hydrothermally from gel precursors in an externally heated Morey-type pressure vessel with an internal silver tubing [31,81].Synthesis conditions were 400 • C at 1 kbar P(H 2 O) during four weeks.The sample was then Na + -saturated at room temperature by contact with aqueous solutions of NaCl (1 mol/L) as described by Dazas et al. [31] Oriented slide was prepared by drying at room temperature an aqueous hectorite dispersion on a glass slide.Hectorite sample was initially equilibrated at 97% relative humidity (RH) over a saturated CuSO 4 solution before being transferred to the chamber and XRD patterns, leading to almost homogeneous 2W and 1W states, were collected at 80% and 28% RH, respectively [31].At 80% RH, the dominant layer is the 2W state (97%) against 3 % for the 1W and 0% for the 0W.At 28% RH, the dominant layer is the 1W state (97%) against 3 % for the 2W and 0% for the 0W.
The algorithms developed initially by Sakharov and co-workers were used to calculate theoretical XRD profiles over the 2-50 • 2θ CuKα range with a trial-and-error approach [82,83].Instrumental and experimental factors such as horizontal and vertical beam divergences, goniometer radius, length, and thickness of the oriented slides were measured and introduced without further adjustment.The mass absorption coefficient (µ*) was set to 45 cm 2 •g −1 , as recommended [84].Additional parameters include the layer-to-layer distance set as proposed by Dazas et al. [31] and the coherent scattering domain size along the c* axis, which was characterized by a maximum value, set to 50 layers, and by a variable mean value (N) [85].N was taken equal to 6 and 5 in the 2W and the 1W state respectively.The z-coordinates of all atoms building up the 2:1 layer were set as determined for hectorite [86].
Calculations of XRD profiles were performed taking into account the hydration heterogeneities (proportion of 0W, 1W, and 2W layers) reported above.The interlayer species density profiles calculated from MD simulations were introduced for the dominant layer according to the methodology proposed elsewhere [7,39,40], whereas the interlayer configurations of the other layer types were set as described by Dazas et al. [31].

Dry Hectorites
The unit cell parameters obtained for the dry hectorite are summarized in Table 5.Few experimental data are available for hydroxylated hectorites containing other cations than Na + .However, the dimensions of the unit cells and the layer-to-layer distances can be compared with talc and other smectites.a, b and γ parameters which correspond to the dimensions of the unit cell in the plane of the clay sheet compare well with the experimental values obtained on talc.The c, α and β parameters are not directly comparable since their value will depend strongly on the arrangement of the counter-ions in between the clay layers.However, the layer-to-layer distance h obtained for Na-hectorite is identical to the value obtained by Bowers et al. [64] and Dazas et al. [31] from XRD patterns, h = 9.7 Å.The layer-to-layer distances for the other hectorites are close to the experimental values obtained on swelling clays with similar charges: 10.7-11.5 Å for smectites containing Cs + [60,61,88], 9.55-10.0Å for smectites containing Ca 2+ [40,62,63], 9.8-10 Å for smectites containing Sr 2+ [40,62,63].
As for talc, the hexagonal cavity at the clay surface obtained with PIM is deformed as seen on Figure 1.As discussed elsewhere [36], the shape of the cavity could result in a smaller space for the cations, which may enter less deeply in the cavity.It explains the larger values of h found with PIM than with clayFF for the small cations (9.6 Å for Na and Ca-hectorite [48]), which cannot enter the cavity in the PIM case.

Hydrated Hectorites
Atomic density profiles along Z-axis obtained using PIM and ClayFF force field and made symmetric with respect to the mid-plane are compared in Figure 3.In the case of 1W hydration state, the interlayer distribution of water is similar, but the equilibrated positions for sodium cations display contrasted behavior.Using the PIM force field, interlayer cations are indeed found to remain on the interlayer mid-plane while they tend to be located closer to the clay surface for ClayFF model.For the 2W hydration state, the situation differs as sodium profiles are similar in both models.The interlayer equilibrium positions for water, however, are significantly different and display an increased broadening for the PIM model compared to the clayFF force field, which is clearly seen on the electronic density profiles of Figure 5.Let us note that the density profiles obtained with ClayFF are every similar to the ones obtained in a previous study with ClayFF and the flexible SPC model [10].
Assessment of the obtained MD configurations from the comparison between experimental and MD-generated XRD profiles is shown in Figure 4 using the structure parameters detailed in the section "XRD materials and methods".For this purpose, the interlayer atomic density profiles from MD simulations are divided into a series of individual atomic planes (separated by a δZ of ∼0.05 Å) containing O, H, or Na as neutral atoms.No thermal fluctuation parameter is considered for these planes because the positional disorder related to temperature is accounted for in the optimized atomic configurations extracted from MD simulations.Comparison between MD-generated and experimental XRD patterns for the 1W hydration state show a good agreement for both PIM and ClayFF models despite the aforementioned difference in atomic Z-positions for cations (Figure 4).Because of the interaction of X-rays with electrons, XRD is extremely sensitive to the overall electronic density distribution and not to the nature of atomic species.The electronic density profile related to both models and compared in Figure 5 is calculated by summing, for each interlayer plane at a given Z-coordinate, the electronic content related to the amount and nature of atoms lying on this plane.As shown in Figure 5, the electronic density profile is principally governed by the distribution of O atoms and even though Na atoms contain a higher amount of electrons than O (i.e., 11 and 8 electrons for Na and O, respectively), the low Na content does not imply significant changes on the electronic profile, which are similar for both PIM and ClayFF models.More differences between MD-generated XRD profiles using PIM or ClayFF models can be noticed in the case of the 2W hydration state, in particular for the intensity ratio between the 007 and 008 reflections (Figure 4).Whereas a good fit is received for all other 00l reflections, the 008:007 intensity ratio is larger for PIM than for ClayFF model.As discussed previously by Ferrage et al. [7], the intensity ratio between these two high-order reflections strongly depends on the heterogeneous distribution of water molecules along Z-axis, an increase in positional disorder inducing a rise on the 008:007 intensity ratio.In agreement with this statement, the comparison of electronic density profile (Figure 5) for the 2W state shows that PIM model induces a more disordered organization of water molecules compared to the ClayFF force field, likely responsible for the slightly better reproduction of the experimental feature with this model.

Beyond Validation: Dynamics
To go further, in this paragraph we compare the dynamical properties obtained with clayFF and PIM.
The diffusion coefficients D were calculated as the slope of the parallel mean-squared displacement of the mobile species ( x 2 + y 2 ) as a function of time in a region where it is linear.Because of the low number of cations in the simulation box, the region was limited to  ps for Na + .For water molecules, the region between 100 and 300 ps was chosen.
We also calculated the intermittent residence time of water molecules around sodium cations.It consists in calculating the survival probability for a molecule to be found within a given distance d of the Na + : P(t) = S I (t) S I (0) where S I (t) = 1 if the water molecule is present in the sphere of radius d centered around the sodium both at time 0 and t, and S I (t) = 0 otherwise.Such a definition allows for molecules to leave and come back into the sphere between 0 and t and characterizes an intermittent survival probability.d was taken as the position of the first minimum of OW-Na radial distribution functions, i.e., 3.3 Å.As a first approximation, P(t) can be fitted by a decreasing exponential exp(−t/τ) where τ can be assimilated to a residence time within the sphere.The water diffusion coefficients were found to be 1.2 × 10 −9 m 2 •s −1 and 9.4 × 10 −10 m 2 •s −1 for clayFF and PIM respectively: water diffuses around 20% more slowly with PIM than with clayFF.This result is encouraging since we found previously that in swelling clays, simulated water diffusion coefficients could be too high when compared to quasi-elastic neutron scattering data [9,32,89].The sodium diffusion coefficients were found to be very close with the two force fields (5 × 10 −10 and 4.5 × 10 −10 m 2 •s −1 for clayFF and PIM respectively).A previous simulation of Na-hectorite with a close water content (10 H 2 O per cation) gave smaller diffusion coefficients, with clayFF and flexible SPC water model: 1.5 × 10 −10 m 2 •s −1 for Na + and 4.8 × 10 −10 m 2 •s −1 for water [10].These differences could come from the water model but also from the layer-to-layer distance, which was 15.1 Å in this study, around 0.4 Å smaller than in this work.Indeed it has been already shown that the layer-to-layer distance could greatly influence the diffusion coefficients [32].Surprisingly, the intermittent residence time of water molecules around sodium cations was found to be smaller with PIM (48 ps) than with clayFF (80 ps).It can be concluded that the slowing down of

Figure 1 .
Figure 1.On the left: view of the hectorite sheets as simulated.h is the layer-to-layer distance.On the right: View of the surface of an hectorite sheet simulated with Polarizable Ion Model (PIM).Yellow: Si; Red: O; White: H; Pink: Mg; light blue: Li.; Dark blue: counter-ion

Figure 2 .
Figure 2. Forces for each atom for one of the hydrated hectorite configurations.The predictions of the classical force field (black lines) for the force (F x , F y and F z ) and dipole (µ x , µ y and µ z ) components are compared to the DFT results (red lines).The atom indices are grouped by atoms types: cations (Na, Mg, Li, Si), clay oxygen atom (O), clay oxygen atoms in hydroxyl groups (Oh) and water oxygen atoms (OW).The groups are delimited by the dashed lines.

Figure 3 .
Figure 3. Atomic density profiles (arbitrary units) of oxygen (red), hydrogen (gray), and sodium (pink) atoms generated from MD simulations using PIM (left) or ClayFF (right) models for monohydrated (1W-top) and bihydrated (2W-bottom) hydration states.Atomic Z coordinates (in Å) are given relative to the mid-plane of the interlayer.

Figure 4 .
Figure 4. Experimental (black crosses) and calculated (red lines) intensities of 00l reflections for hectorite under monohydrated (1W-left) and bihydrated (2W-right) state and using PIM (top) or ClayFF (bottom) models.The vertical gray bars indicate a modified scale factor for the high-angle regions as compared to the low-angle part of the patterns.00l reflections are indexed on the top part of the figure.For 2W state, a zoom of the 007 and 008 reflections is shown for both models.

Figure 5 .
Figure 5. Interlayer electronic density profiles (arbitrary units) calculated from atomic density profiles generated using PIM (gray) or ClayFF (black) models for monohydrated (1W-left) and bihydrated (2W-right) hydration states.Atomic Z coordinates (in Å) are given relative to the mid-plane of the interlayer.
2−, OH − , Mg 2+ , Al3+and Si 4+ contrary to non polarizable force-fields where partial charges are used to describe charge transfers in ionocovalent systems.Only the charge transfer within the hydroxyl group, of total charge −1, is modelled by partial charges on the corresponding atoms: O

Table 1 .
Characteristics of the Dang-Chang water model.The oxygen atoms interact via a Lennard-Jones potential.The negative charge and the dipole are carried by the fourth site MW.

Table 2 .
χ 2 for the dipoles and the forces for dehydrated and bihydrated hectorite.

Table 3 .
Parameters of the PIM force field for the interactions between lithium and the other atoms.

Table 4 .
Characteristics of the simulations boxes for hydrated hectorites.

Table 5 .
Unit cell parameters of dry hectorites with the PIM model.The error on the last significant digit is reported between parenthesis.