Temperature dependence study of water dynamics in Fluorohectorite clays using Molecular dynamics simulations

In this work, we have carried out molecular dynamics (MD) simulation techniques to study the diffusion coefficient of interlayer molecules at different temperature. Within the wider context of water dynamics in soils, and with a particular emphasis on clays, we present here the translational dynamics of water in clays, in a bi-hydrated states. We focus on temperatures between 293 K and 350 K, i.e., the range relevant to the environmental waste packages. A natural hectorite clay of interest is modified as a synthetic clay, which allows us to understand the determinantal parameters from MD simulations through a comparison with the experimental values. The activation energy E a determined by our simulation is [8.50 - 16.62] KJ/mol. The calculated diffusive constants are in the order of 10^{-5} {cm^{2}}/s. The simulation results are in good agreement with experiments for the relevant set of conditions, and they give more insight into the origin of the observed dynamics.


Introduction
Clays are characterized by high-surface area, thus, making them a promising material for applications involving the adsorption of other molecules on the surface.Because of the layered nature of these materials, adding various molecules between layers might result in adsorption with significant binding energies.Additionally, clays are plentiful, nontoxic, inexpensive, and reusable [1].
When clay minerals are used, two important concerns develop.These are the depletion of natural deposits, particularly, those with easy mining access, and the presence of natural clays in many phases rather than in a pure single phase.To solve these difficulties, in recent decades, geologists, material scientists, and geochemists have developed a great interest in the laboratory synthesis of clays [11][12][13].In this study, we make a theoretical analysis of the properties of dynamics of water in fluorohectorite, a synthetic smectite group of clay materials.Fluorohectorite is considered to be an important material to the clay science, supporting a claim by scholars which stated that "Clays may be considered as the material of the 21 st century" [14].Analyzing water dynamics in anisotropic and confined clay minerals could help to understand transport properties in related materials such as split solids, porous media, and soil models.Water is important for clay materials when it comes to the storage of radioactive wastes since the waste is surrounded by water and ions that can spread inside the clay.As a result, an understanding of the behavior of water dynamics on various length scales as well as temperatures is necessary.In this work, we conduct computer simulation to explore the water dynamics in fluorohectorite clay interlayers at temperature ranges of 293-350 K.The paper is organized as follows.In the next section (Section 2), a detail account of the computational method is presented.Results and discussion are presented in Section 3, with the conclusion being presented in Section 4.

Computational Methods
where M denotes some form of an intercalated cation between clay layers.In the present work, the intercalated cations considered are monovalent cations such as Li þ , Na þ , K þ , and Cs þ .The simulation box is shown in Figure 1, and the corresponding crystallographic positions of the atoms in the clay unit cell are given elsewhere [15].Our periodic simulation cell has three dimensional periodic boundary conditions.The sheets of clay are considered as rigid molecules.The van der Waals interactions between water-water, water-clay, watercation, clay-clay, clay-cation, and cation-cation are defined by the Lennard-Jones (LJ) potential model.Summations across all interaction sites generate the system's total potential energy, where the pair-wise interaction is given by where q i and q j are the partial charges carried by the atoms; and ϵ ij and σ ij are the LJ parameters obtained from the corresponding atomic parameters using the Lorentz-Berthelot mixing rule [16]: and All of these parameters are determined by the type of force field used to characterize the system.In our investigation, we chose a clayFF force field that is suitable to our synthetic clay since it has been found to provide transport properties of water, that are close to experiment (see elsewhere in Section 3).The water molecule in this clay-water interaction is represented by the SPC water model.This model is favored over the TIP4P/2005 model for H 2 O molecules [17].The distribution of interlayer cation and water perpendicular to the plane of the clay layers are analyzed using atomic density profile.It is used to investigate the extent of presence of the various species in the space between the layers of clay.The local coordination distributions between interlayer ion and the basal oxygen atoms of the clay layer (OB) or water molecules (OW, HW) are characterized by radial distribution function (RDF).The RDF for species Y around species X is calculated as follows: where ρ Y is a density of species Y, the fraction dðN X−Y Þ dr is the average number density of a particle of species Y lying in the region r to r þ dr around a particle of species X, and N X−Y is the coordination number for species Y around species X.Furthermore, the RDF provides the probability of finding a particle at a certain distance from a reference particle.But, there is no probability of finding an atom from the first neighbor shell within up to the distance of the atom diameters.This region of zero RDF is known as exclusion region [18].The mean square displacement (MSD) of the ion and water molecule during the production run is analyzed to determine the dynamical characteristics of the interlayer cation and water.To conduct the molecular dynamics simulations as well as make analysis of quantities such as the self diffusion coefficients, the GROMACS package is utilized.The executable "gmx msd" generates an MSD data file about average MSD as a function of time.Then we acquired the slopes of the MSD graphs by doing a linear fit on this data.Finally, as formulated by Einstein relation, dividing the slope by the factor 6, see Equation ( 6), we get self diffusion coefficients of interlayer cations and water molecules as follows where N m is the number of selected species, r j ðtÞ is the center of mass position of the j th species at time t.The horizontal/ lateral diffusion coefficient is calculated using, where the slope of the MSD parallel to the clay (xy plane) is a function of the time considered (up to 1,000 ps in this study).Some literature make a direct link between the diffusion coefficient and temperature, given by Arrhenius equation [19], as follows: ) The density profile of interlayer cation and water molecules of the different fluorohectorite clays (MFht; M = Li, Na, K, Cs) as a function of z coordinate at different values of temperature.From left to right, the temperatures are 293, 300, 323, and 350 K, respectively.The top side is for the ions (M) while the bottom side is for water.Colors: blue color and solid line-LiFht, orange color and dashdot-NaFht, green color and dashed line-KFht, and black color and solid line-CsFht.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.
Advances in Condensed Matter Physics where D 0 is a diffusion coefficient at T → 1, E a is activation energy for diffusion, k B is the Boltzmann constant, and N A is Avogadro number.This equation can be further written as follows: The activation energy, E a , can be obtained by a linear fit of ln ðD=D 0 Þ vs. 1  T according to Equation (9), by taking the slope of this graph and then multiplied by ð − N A k B Þ.

Molecular Dynamics Simulation.
The steepest descent algorithm is used upon optimizing the structures and energies.The system with initial atomic positions assigned within the unit cell is optimized until the net force on the atoms falls below the threshold value and the total energy is converged within tolerance value.The energy minimization process is repeated twice while maintaining other circumstances and parameters the same for analyzing the diffusion phenomena.When a flexible bond is used to allow atoms to move apart from one another in a controlled way, a restricted bond is used to ensure that the new constrained positions do not experience strong forces.After energy minimization, the system is brought into temperature and pressure equilibration.The system's thermodynamic characteristics change depending on many factors such as temperature, pressure, density, etc.The accuracy of thermodynamic properties to be determined would be impacted by the changes to these factors.Consequently, a system must undergo equilibration before a production run [20].The clay-water system is equilibrated at four different temperatures ranging from 293 to 350 K, and isobaric pressure of 1 bar, with isothermal compressibility of 4.5 × 10 −5 /bar.For a temperature coupling, a velocity rescale thermostat is employed, and for pressure coupling, a Berendsen barostat is used.The coupling time constants for the thermostat and the barostat are 0.1 and 2.0 ps, respectively.The duration of the equilibration is 1,000 ps.
The Particle Mesh Ewald (PME) summation algorithm is used for long range interactions.The cutoff parameter of 1.2 nm is taken with periodic boundary conditions for Coulomb and LJ interactions.After equilibration, production run is carried out to calculate the thermodynamic properties of the system such as partial density (ρðrÞ), RDF (gðrÞ), and diffusion coefficient (D), using NVE ensemble.All structural and dynamical quantities are performed for 1,000 ps with the time step of 0.001 ps.

Results and Discussion
The presentations in this section are a revised version of the outcomes reported in preprint, which are available in arXiv repository [21].For the analysis of the results, our approach resembles the ones adopted in a literature [22].Figure 2 shows the calculated atomic density profiles of the interlayer ions and water at different temperature inside the pores of LiFht, NaFht, KFht, and CsFht clays.The profiles were computed along a direction perpendicular to the clay surface and averaged over the interlayer region of the simulation box using a 1,000 ps production time frame.A detailed analysis requires the RDF of interlayer cation ions with respect to reference oxygen atoms of water (OW) are shown in Figures 3-6, with red color.For K and Cs-Fluorohectorite, the 1 st neighboring peak has gðrÞ peak intensity values of 7.64 at T = 293 K; 7.30 at T = 300 K; 7.22 at T = 323 K; 7.33 at T = 350 K; 5.47 at T = 300 K; 5,46 at T = 300 K; 5.45 at T = 323 K; and 5.38 at T = 350 K.However, for Li and Na-Fluorohectorite, the 1 st neighboring peak has gðrÞ peak intensity values of 25.13 at T = 293 K; 24.50 at T = 300 K; 23.92 at T = 323 K; 22.45 at T = 350 K; 19.76 at T = 293 K, 19.55 at T = 300 K; 18.47 at T = 323 K; and 17.24 at T = 350 K.This result simply shows that the peak intensity change owing to the presence of different interlayer cations.In addition to this, comparing with the hydration shell around K þ and Cs þ cations, we found that the 1 st neighbor cation-oxygen are assembled more densely than those in the and Cs þ , respectively.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.4 Advances in Condensed Matter Physics case of Li þ and Na þ cations.The RDF of interlayer cation with respect to basal oxygen of clay (OB) and with respect to hydrogen of water (HW) are represented by black and blue colors, respectively, as shown in the Figures 3-6.We can see the exclusion region up to a certain distance, where the probability of finding the particle with respect to the reference particle is zero.At a large value of r, the correlation between them decreases and eventually there would not be any long range correlation, i.e., gðrÞ → 1 [23].The exclusion regions for cation-OW and cation-OB correlation pairs are provided in Table 1 at all temperature.The coordination number of interlayer cations is the number of molecule of water and basal oxygen of the clay in the immediate neighbor of the cations.It depends on the distance between the water molecule as well as the basal oxygen of clay and cations.In Tables 2 and 3, the coordination number of each correlation pair is presented.
Figure 7 shows the MSD plots of interlayer cations (Li þ , Na þ , Cs þ , and K þ ), and H 2 O molecule in LiFht, NaFht, KFht, and CsFht, at different temperatures ranging from 293 to 350 K, respectively.The corresponding simulated results of self-diffusion coefficients in 3D and 2D are shown in Tables 4 and 5, respectively.It is clearly observed that when temperature increases, the generated velocities of the ions and water molecules also increase and the density of the system decreases.This provides more space for the molecule to execute random walk.Due to this, the MSD increases.From Einstein's relation in Equation ( 6), as the MSD increases, the self-diffusion coefficient also increases.
In the confining environment of nanochannels, the lateral self diffusion coefficients for motion parallel to the clay sheets are more relevant and necessary.Due to this, we have calculated the lateral diffusion coefficient of both interlayer cation and water for all temperature values and the different types of fluorohectorite clay.The result is presented in Table 5.According to a literature [25], the lateral diffusion coefficient of water in Na-montmorilonite is D xy ¼ ½7:7 AE 1:0 × 10 −6 cm 2 s −1 .Generally, in all cases, the self-diffusion coefficient of both cations and water depend on temperature.As the temperature increases, the diffusion coefficients for water seem to increase.Such trend is also reported in a literature [22], with the reported values of the diffusion coefficients of water being closely comparable to the investigations in this work.Figure 7 shows the of MSD of ions and water at different temperatures.It looks that the MSD and the diffusion coefficients of water is greater than that of the cations.The main factor which affects the diffusion coefficient of interlayer water and cations seems to be the electrostatic interaction between water and cations of the clay sheet (OW, HW).Experiment literature [5] indicated formation of hydrogeneous structure with cation water (M-H 2 O) complex in the interlayer regions.In support of this, this work has investigated that a pathway for the mobility and retention of water in the clay interlayers is through a bond formation and bond breaking in a cation-water complex (analyzed as OW, HW).
Figures 8 and 9 shows the Arrhenius plot of the diffusion coefficient plot of interlayer molecules.The pre-exponential factor (D 0 ) can be determined by extrapolating the graph to zero.The activation energy for the diffusion of water molecule as calculated from the simulation is 0.14, 0.13, 0.09, and 0.12 eV per molecule in LiFht, NaFht, KFht, and CsFht, respectively, and the corresponding activation energy for the diffusion of Li þ , Na þ , K þ , and Cs þ are 0.17, 0.10, 0.11, and 0.09 eV per ion, respectively, in 3D motions.

Conclusion
The transport characteristics of cations and water molecules in bihydrated LiFht, NaFht, KFht, and CsFht clays shows that Na þ and Li þ ions exhibit significant diffusion motion compared to K þ and Cs þ ions, while the latter exhibits a more constrained motion as well as better potential for screening negative charges (wastes).Water seems to have relatively larger diffusion coefficients in KFht and CsFht clays when compared to the values in LiFht and NaFht clays.Such a higher diffusion prospect of water may mean a more suitable prospect for the applications in wastewater treatment.The self-diffusion coefficients calculations reveal that the values increase as the temperature increases.Diffusion coefficients for bilayer states are typically in the order of 10 −5 cm 2 s −1 in most neutron scattering experiments on natural clays, like montmorillonite, hectorite, or vermiculite.In our synthetic hectorite clay, called fluorohectorite clays, the outcomes from NaFht, LiFht, KFht, and CsFht also appears to be on the order of 10 −5 cm 2 s −1 .Such a comparable diffusion coefficients in the different fluorohectorite clays may also mean competitive applications for the absorption of wastewater.A lateral (2D) dynamics activities contribute significantly to the observed 3D properties, which may mean that the clay interlayer is a crucial pathway for the transportation of ions and water.As a future extended activity to the present study, one may consider the possibility to analyze data by including several effects such as rotation and vibrations; as well as by using different models for water.

FIGURE 3 :
FIGURE 3:  The RDF analysis of gðrÞ at T = 293 K. Colors: red for cation-OW, black for cation-OB, and blue for cation-HW.Cations from left to right are Li þ , Na þ , K þ , and Cs þ , respectively.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.

25 FIGURE 5 :
FIGURE 5:  The RDF analysis of gðrÞ at T = 323 K. Colors: red for cation-OW, black for cation-OB, and blue for cation-HW.Cations from left to right are Li þ , Na þ , K þ , and Cs þ , respectively.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.

FIGURE 4 : 15 FIGURE 6 :
FIGURE 4:  The RDF analysis of gðrÞ at T = 300 K. Colors: red for cation-OW, black for cation-OB, and blue for cation-HW.Cations from left to right are Li þ , Na þ , K þ , and Cs þ , respectively.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.

FIGURE 7 :
FIGURE 7:  The mean square displacement graph as a function of time of interlayer water.The counterpart for the interlayer cations is shown in insets.Colors: blue color and solid line for T = 293 K, orange color and dashed line for T = 300 K, green color and dotted line for T = 323 K and black color and dashdot line for T = 350 K.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.

FIGURE 8 :FIGURE 9 :
FIGURE 8: The (a) 3D coefficient vs. temperature graph and (b) The natural logarithm of 3D diffusion coefficient vs. inverse temperature graph of interlayer ions of LiFht, NaFht, KFht and CsFht clays.(c) The natural logarithm of 3D diffusion coefficient vs. inverse temperature graph of interlayer water molecule.Colors: blue-LiFht, orange-NaFht, green-KFht, and black-CsFht.For the interpretations of the references to color in this plot legend, the reader is referred to the web-version.

TABLE 1 :
The exclusion region of each correlation pair in Å for different temperatures.

TABLE 2 :
Coordination numbers in LiFht and NaFht clays at different temperature.

TABLE 3 :
Coordination numbers in KFht and CsFht clays at different temperature.