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

Abstract

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 field 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 field are in very good agreement with experimental results. In the bihydrated state, the less structured electronic density profile obtained with PIM compared to the one from the state-of-the-art non-polarizable force field clayFF explains the slightly better agreement between the PIM results and experiments.

1. Introduction

Clays are widely used in industrial applications and environmental contexts because of their remarkable retention properties. Molecular simulations provide a detailed understanding of the system at the atomic scale and constitute a complementary tool to experiments [1,2,3,4,5,6,7,8,9,10,11,12], analytical theories and modeling at larger scales [11,13,14,15,16,17,18,19,20,21]. They are based on the calculation of the interatomic interactions which are defined by a set of equations and parameters, the force-field. Several force-fields have been developed for clay/water systems which offered a realistic description of such systems at the atomic scale [22,23,24,25,26,27].

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.

2. Parametrization of PIM Force Field

2.1. Characteristics of the Polarizable Ion Model

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

$${\mathrm{V}}_{\mathrm{total}}={\mathrm{V}}_{\mathrm{Charge}}+{\mathrm{V}}_{\mathrm{Dispersion}}+{\mathrm{V}}_{\mathrm{Repulsion}}+{\mathrm{V}}_{\mathrm{Polarization}}.$$

(1)

The charge term corresponds to the Coulomb interaction (here in atomic units) between two atoms,

$${\mathrm{V}}_{\mathrm{Charge}}=\sum _{i<j}\left(\frac{q_i{q}_j}{r_{ij}}\right),$$

(2)

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${}^{2-}$, OH${}^{-}$, Mg${}^{2+}$, Al${}^{3+}$ 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${}_{\mathrm{OH}}^{(2-\delta )-}$ and H${}_{\mathrm{OH}}^{(1-\delta )+}$ with $\delta =+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]:

$${\mathrm{V}}_{\mathrm{Dispersion}}=-\sum _{i<j}\left[f_6^{ij}\left(r_{ij}\right)\frac{C_6^{ij}}{{\left(r_{ij}\right)}^6}+f_8^{ij}\left(r_{ij}\right)\frac{C_8^{ij}}{{\left(r_{ij}\right)}^8}\right],$$

(3)

where $C_6^{ij}$ and $C_8^{ij}$ are the dipole-dipole and dipole-quadrupole dispersion coefficients. The Tang-Toennies damping function $f_n^{ij}$ is used to correct the short-range interaction as [53]:

$$f_n^{ij}\left(r_{ij}\right)=1-e^{-b_n^{ij}{r}_{ij}}\sum _{k=0}^n\frac{{(b_n^{ij}{r}_{ij})}^k}{k!},$$

(4)

where $1/b_n^{ij}$ is the range of the damping.

The repulsion term in Equation (1) is modelled as:

$${\mathrm{V}}_{\mathrm{Repulsion}}=\sum _{i<j}{A}_{ij}{e}^{-B_{ij}{r}_{ij}},$$

(5)

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:

$${\mathrm{V}}_{\mathrm{Polarization}}=\sum _{i<j}\left[\frac{q_i{\mathit{r}}_{ij}·{\mathbf{\mu }}_j}{{r_{ij}}^3}{g}_4^{ij}\left(r_{ij}\right)-\frac{{\mathbf{\mu }}_i·{\mathit{r}}_{ij}{q}_j}{r_{ij}^3}{g}_4^{ji}\left(r_{ij}\right)+\frac{{\mathbf{\mu }}_i·{\mathbf{\mu }}_j}{r_{ij}^3}-\frac{3\left({\mathit{r}}_{ij}·{\mathbf{\mu }}_i\right)\left({\mathit{r}}_{ij}·{\mathbf{\mu }}_j\right)}{r_{ij}^5}\right]+\sum _i\frac{{\left|{\mathbf{\mu }}_i\right|}^2}{2{\alpha }^i},$$

(6)

where ${\alpha }^i$ is the polarizability of ion i, ${\mathbf{\mu }}_i$ and ${\mathbf{\mu }}_j$ are the induced dipoles, $g^{ij}$ is the short-range correction to the multipolar expansion by the Tang-Toennies damping function:

$$g_4^{ij}\left(r_{ij}\right)=1-c_{ij}{e}^{-b_D^{ij}{r}_{ij}}\sum _{k=0}^4\frac{{(b_D^{ij}{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 dynamic step by minimizing the polarization energy:

$$\left(\frac{\partial {\mathrm{V}}_{\mathrm{Polarization}}}{\partial {\mu }_{\alpha }^i}\right)=0.$$

(8)

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_6^{ij}$, $C_8^{ij}$ and $b_n^{ij}$ 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.

2.2. Optimization Procedure

The optimization procedure aims at finding the set of parameters ($A_{ij}$, $B_{ij}$, $c_{ij}$ and $b_D^{ij}$) 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_D^{ij}$ 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:

• Generation of a series of representative configurations using classical molecular dynamics (MD)
• 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 ${\chi }_{\mathrm{Dipole}}^2$ with respect to the parameters of the polarization term (V${}_{\mathrm{Polarization}}$) and of ${\chi }_{\mathrm{Force}}^2$ with respect to the repulsion term (V${}_{\mathrm{Repulsion}}$):

  $${\chi }_{\mathrm{Dipoles}}^2(b_D^{ij},c_{ij})=\frac{1}{N_{\mathrm{conf}}}\frac{1}{N_{\mathrm{atom}}}{\sum }_{\mathrm{conf}}{\sum }_{\mathrm{atom}}\frac{\left|\right|{\mathbf{\mu }}^{\mathrm{classical}}-{\mathbf{\mu }}^{\mathrm{DFT}}{\left|\right|}^2}{\left|\right|{\mathbf{\mu }}^{\mathrm{DFT}}{\left|\right|}^2},$$

  (9)

  $${\chi }_{\mathrm{Forces}}^2(A_{ij},B_{ij},b_D^{ij},c_{ij})=\frac{1}{N_{\mathrm{conf}}}\frac{1}{N_{\mathrm{atom}}}{\sum }_{\mathrm{conf}}{\sum }_{\mathrm{atom}}\frac{\left|\right|{\mathit{F}}^{\mathrm{classical}}-{\mathit{F}}^{\mathrm{DFT}}{\left|\right|}^2}{\left|\right|{\mathit{F}}^{\mathrm{DFT}}{\left|\right|}^2},$$

  (10)

  where $N_{\mathrm{conf}}$ is the number of configurations on which DFT calculations are performed, $N_{\mathrm{atom}}$ is the number of atoms per configuration, ${\mathbf{\mu }}^{\mathrm{classical}}$ are the dipoles obtained by classical molecular dynamics using a given set of parameters. ${\mathit{F}}^{\mathrm{classical}}$ are the forces obtained by classical molecular dynamics.

2.3. Simulation Details

2.3.1. 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_{\mathrm{Li}-\mathrm{Cation}}$ and $b_D^{\mathrm{Li}-\mathrm{Cation}}$, 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].

2.3.2. 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. 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.

${\mathit{d}}_{\mathbf{OH}}$ (Å)	${\mathit{d}}_{\mathbf{OM}}$ (Å)	Angle (${}^{\circ }$)	${\mathit{ϵ}}_{\mathbf{O}}$ (kcal/mol)	${\mathit{\sigma }}_{\mathbf{O}}$ (Å)	${\mathit{q}}_{\mathbf{H}}$	${\mathit{\alpha }}_{\mathbf{M}}$ (Å${}^3$)
0.9752	0.215	104.52	0.1825	3.2340	0.5190	1.444

.

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].

2.4. 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, ${\chi }_{\mathrm{Dipoles}}^2$ and ${\chi }_{\mathrm{Forces}}^2$ are shown in

Table 2. ${\chi }^2$ for the dipoles and the forces for dehydrated and bihydrated hectorite.

Systems	${\mathit{\chi }}_{\mathbf{Dipoles}}^2$	${\mathit{\chi }}_{\mathbf{Forces}}^2$
Dry hectorite	0.123	0.561
Wet hectorite	0.170	3.025

. 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. ${\chi }_{\mathrm{Dipoles}}^2$ and ${\chi }_{\mathrm{Forces}}^2$ 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. ${\chi }_{\mathrm{Dipoles}}^2$ and ${\chi }_{\mathrm{Forces}}^2$ 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. Parameters of the PIM force field for the interactions between lithium and the other atoms. The parameters $C_6^{ij}$, $C_8^{ij}$ and $b_n^{ij}$ between Li and O are equal to the interactions O${}_{\mathrm{OH}}$-X (with X = Al, Mg, Si, O or O${}_{\mathrm{OH}}$) [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.

Damping Interaction between ${\mathit{q}}_{\mathit{j}}$ and ${\mathit{\mu }}_{\mathit{i}}$
Ion Pair (ij)	${\mathit{A}}_{\mathit{ij}}$ (Ha)	${\mathit{B}}_{\mathit{ij}}$ (Å−1)	${\mathit{C}}_{\mathbf{6}}^{\mathit{ij}}$ (Ha·Å6)	${\mathit{C}}_{\mathbf{8}}^{\mathit{ij}}$ (Ha·Å8)	${\mathit{b}}_{\mathit{n}}^{\mathit{ij}}$ (Å−1)	${\mathit{c}}_{\mathit{ij}}(-)$	${\mathit{b}}_{\mathit{D}}^{\mathit{ij}}$ (Å−1)
Li-O	9.68	2.78	0.0477	0.156	4.17	1.19	2.79
Li-O${}_{\mathrm{OH}}$	24.6	3.82	0.0477	0.156	4.17	1.16	2.78
Li-OW	0.427	3.50	0.0477	0.156	4.17	-	-
Li-MW	-	-	-	-	-	3.75	1.24
Li-Na	-	-	-	-	-	0.587	4.64
Li-Cs	-	-	-	-	-	0.0423	4.99
Li-Ca	-	-	-	-	-	3.87	3.64
Li-Sr	-	-	-	-	-	0.895	4.69

.

3. 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.

3.1. 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. Characteristics of the simulations boxes for hydrated hectorites.

System	Super Cell Dimensions	A (Å)	B (Å)	h (Å)	Number of H${}_2$O per Unit Cell
monohydrated (1W)	5 × 2 × 2	26.3	18.2	12.53	4.2
bihydrated (2W)	5 × 2 × 2	26.3	18.2	15.51	8.9

.

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.

3.2. 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 ${}^{\circ }$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${}^{\circ }2\theta$ CuK$\alpha$ 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 ($\mu$*) 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].

3.3. Results and Discussion

3.3.1. Dry Hectorites

The unit cell parameters obtained for the dry hectorite are summarized in

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

Counterion	a (Å)	b (Å)	c (Å)	h (Å)	$\mathbf{\alpha }$ (deg)	$\mathbf{\beta }$ (deg)	$\mathbf{\gamma }$ (deg)
Na${}^{+}$	5.235 (1)	9.108 (4)	9.813 (7)	9.711 (1)	91.8 (2)	96.4 (2)	90.02 (5)
Cs${}^{+}$	5.258 (5)	9.106 (7)	10.674 (7)	10.524 (5)	89.90 (6)	99.53 (7)	89.95 (4)
Ca${}^{2+}$	5.255 (9)	9.15 (2)	10.23 (1)	9.676 (1)	90.1 (1)	108.1 (2)	90.2 (1)
Sr${}^{2+}$	5.231 (6)	9.11 (2)	10.28 (2)	9.767 (5)	90.1 (2)	106.9 (2)	90.1 (1)
talc (exp) [87]	5.293	9.179	9.469	9.381	90.57	98.91	90.03

. 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 $\gamma$ 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, $\alpha$ and $\beta$ 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.

3.3.2. 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 $\delta 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.

3.4. 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 ($\langle x^2+y^2\rangle$) 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 [50–150] 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\left(t\right)=\frac{\langle S_I\left(t\right)\rangle }{\langle S_I\left(0\right)\rangle }$$

where $S_I\left(t\right)=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\left(t\right)=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\left(t\right)$ can be fitted by a decreasing exponential $exp(-t/\tau )$ where $\tau$ 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 water molecules is less due to the presence of the cations in the case of PIM than in the case of clayFF. The influence of the clay surfaces could be stronger in this case. This behavior could still differ more when changing the nature of the cation (more charged or polarizable cations); this will be the subject of future works.

4. Conclusions

In this work, we extend the PIM polarizable force field to hectorite clays, using the same procedure as previously for pyrophyllite, talc and montmorillonite. The structure of the atomic clay structure in dehydrated, 1W and 2W states was validated by direct comparison with XRD patterns obtained on the same systems. Such direct comparison between experimental and MD-generated XRD patterns shows that the atomic density profiles resulting from our PIM model and the very popular ClayFF model are in overall good agreement with the experimental data, but that the former more accurately reproduces features of the water distribution in the interlayer. Note that this approach is not sufficiently sensitive to assess the relative merits of the two force fields for Na or H atoms, due to their low amount in the unit cell and their low electronic content, respectively. Additional experimental data such as neutron diffraction on deuterated specimen could provide additional constraints on the validation of the atomic distribution of H atoms [7].

We now plan to extend this approach to fluorohectorites. When put in the presence of water, such synthetized clays show much less interstratification than natural samples [31,90]. As a consequence, almost pure mono- and bihydrated states can be obtained from the experimental point of view and direct comparison with molecular dynamics simulations is easier compared to natural clays [89]. The comparison between PIM results and existing structural and dynamical experimental data on this specific clay is the next step towards PIM validation.

We also plan to extend PIM to ions with large radii, the polarizability of which will play a crucial role on sorption properties. As we showed in this paper, the mechanisms of diffusion can be rather different from one force field to the next. The differences between PIM and non polarizable force fields could be enhanced in this case and change our understanding of these systems. This is especially important for applications like soil and ground water remediation or radioactive waste storage.