Effect of Crystal Chemistry Properties on the Distribution Characteristics of H 2 O and Na + in Na-Montmorillonite Interlayer Space: Molecular Dynamics Simulation Study

: At monolayer hydration state, the spatial distribution of H 2 O and Na + in the interlayer of Na-montmorillonite (Na-MT) with different crystal chemistry properties was investigated by the molecular dynamics simulation method. The simulation results show that when layer charge density increases, H 2 O will move and form hydrogen bonds with O in tetrahedral surfaces (O t ) at a distance of 1.676 ± 0.043 Å. The impact of isomorphic substitution on the relative concentration of H 2 O depends largely on the layer charge density of Na-MT, when layer charge density is high, H 2 O move obviously to both sides of Na-MT sheets with the increase of octahedral substitution ratio. Nevertheless, Na + coordinate with O t at a distance of 2.38 Å, and the effect of isomorphic substitution ratio on the diffusion of Na + is opposite to that of H 2 O. The mobility of both H 2 O and Na + decreases with the increase of layer charge density or tetrahedral substitution ratio. The radial distribution function of Na-O w (O in H 2 O) shows that the coordination strength between Na + and O w decreases with the increase of layer charge density or tetrahedral substitution ratio, and Na + are hydrated by four H 2 O at a Van der Waals radius of 2.386 ± 0.004 Å. The research results can provide a theoretical basis for the efficient application of Na-MT at the molecular and atomic levels.


Introduction
Montmorillonite (MT) is a typical layered aluminosilicate clay mineral and is abundant in mud stone and shale rock [1]. The structure of MT is one Al-O(OH) octahedral sheet sandwiched between two Si-O tetrahedral sheets [2]. Isomorphic substitution in the structure of MT is a requirement, that is Al 3+ are substituted by Mg 2+ in the octahedron and Si 4+ are substituted by Al 3+ in the tetrahedron [3]. As a result of heterovalent isomorphic substitution, the layers of MT are permanently negatively charged. Therefore, some exchangeable cations such as Na + , Ca 2+ in the natural environment can enter the layers of MT to balance the negative charges between MT layers and form various types of MT [4,5].
MT has a wide range of applications in many fields due to its special crystal chemistry properties. It can not only be used for deposition of contaminants and radioactive waste because of its low permeability and good radionuclide retention ability [6,7], but also plays significant roles in the pharmaceutical industry, cosmetic products, and mining operations [8][9][10].
The crystal chemistry properties of MT have a great influence on its application performance and have been extensively studied experimentally by Fourier transform infrared (FTIR) spectroscopy, atomic force microscopy (AFM), X-ray diffraction (XRD), and X-ray photoelectron spectroscopy (XPS) techniques et al. in recent years [11][12][13][14][15]. However, the above methods make it difficult to visualize the crystal chemistry properties of MT from the microscopic perspective. It is well known that molecular simulation has become a popular tool to supplement theoretical research for MT. It can not only provide a detailed insight into the structure of MT at the molecular and atomic levels, but also describe the dynamic behavior of particles in the interlayer domain of MT accurately and explain the interaction mechanism of particles in MT [16][17][18].
Numerous works about the hydration mechanism of MT have been successfully undertaken by the molecular simulation method. Yi et al. investigated the Na-MT-water interface system by molecular dynamics simulation and experimental methods, and verified that there was a hydration shell on the surface of Na-MT with the thickness of approximately 16.3 Å [19], which was composed of six ordered water molecule layers [20]. Titiloye, J.O. et al. studied the structure and dynamics of methane in hydrated K-MT and Na-MT under conditions encountered in the sedimentary basin using computer simulation methods, and found that the K-MT surfaces had a higher affinity for methane than that of Na-form [21]. Voora et al. compared different MT with different interlaminar cations (Li + , Na + , K + , Mg 2+ , and Ca 2+ ) by density functional theory (DFT), and found that the basal spacing correlated in a linear manner with the ionic radius of the interlayer cation [22]. Moreover, other researchers also found that cation with weak interaction force, such as Na + , would be hydrated easily [23][24][25]. Zhou et al. simulated the changes of the interlayer structure in Na-MT under the high temperature and pressure of basin conditions. Their results showed that the monolayer hydrate was more stable than the bilayer hydrate at a burial depth of 7.0 km, and with burial depth increased, the basal spacing of the monolayer and bilayer hydrates changed to varying degrees [26]. In addition, Akinwunmi et al. studied the swelling behavior of Na-MT and found that the swelling pressure appeared rapidly due to the temperature rising [27], however, Peng et al. found that it was the exchangeable cations that provided the dominant driving force for MT to swell through a molecular dynamic study [28].
Regarding the impact of crystal structure on the hydration properties of MT, Qiu et al. studied the effect of charge density on the hydration properties of Ca-MT by combining molecular dynamics simulation and experimental methods. We found that the polarity of Ca-MT sheets was strengthened when the charge density increased, so the interlayer particles were more likely to be drawn to tetrahedral surfaces [29], Shen et al. also drew a similar conclusion [30]. Wungu et al. compared the effect of isomorphic substitution or a non-isomorphic substitution on the geometric structure and electronic properties of MT using DFT calculations. They demonstrated that the isomorphic model possessed insulator-like behavior and was more stable than the conventional model based on the sorption energy [31].
It can be seen from above that the related factors (such as the types of interlayer cations, environmental conditions and crystal structure) affecting the hydration properties of MT have been studied by molecular simulation. However, the above literatures did not fully take into account the crystal structure characteristics of Na-MT in nature, such as the characteristics of layer charge density (from Wyoming to Cheto type) and the difference of isomorphic substitution sites in the structure of Na-MT.
In this work, according to the crystal structure characteristics of Na-MT in nature, 20 types of Na-MT crystal models with different crystal chemistry properties were established. Furthermore, a molecular dynamics simulation was performed to study the influence mechanism of layer charge density and isomorphic substitution sites on the spatial distribution of H2O and Na + in the interlayer domain of Na-MT. The research results are useful for understanding the hydration properties of Na-MT with different crystal chemical properties at the atomic and molecular levels.

Model Establishment
The Na-MT models which contain different layer charge density and isomorphic substitution sites were established firstly. The semi-crystal structure formula of these models was Nax + y(Al2 -xMgx)(Si4 -yAly)O10(OH)2 nH2O, where x + y is the layer charge density, x is the number of magnesium atoms in the octahedron that replace aluminum atoms, y is the number of aluminum atoms in the tetrahedron that replace silicon atoms. The space group of Na-MT is C2/m [32]. The fundamental parameters of the Na-MT model in the simulation were originated from empirical models introduced by Chang et al. [33]: a = 0.520 nm, b = 0.897 nm, the value of c depends on the amount of H2O, α = γ = 90.0°, and β = 99.0°. The Na-MT model consists of two TOT layers, each containing replicated pyrophyllite cells in the form of 4a × 4b × 1c along x, y, and z directions, respectively.
According to the research purposes and crystal structure of natural Na-MT [29], the established semi-cell layer charge of Na-MT was determined as 0.38, 0.50, 0.58 and 0.70 in this paper. The corresponding numbers of Na + which balance the negative potential of Na-MT layers is 6, 8, 9 and 11. The four types of Na-MT were defined as Na-MT1, Na-MT2, Na-MT3 and Na-MT4, respectively. The isomorphic substitution mainly occurs in octahedron (Mg 2+ substitutes Al 3+ ) and tetrahedron (Al 3+ substitutes Si 4+ ), but always obeys Loewenstein's rule that the substitution sites are random and two adjacent positions can not be replaced at the same time [34]. The ratio of octahedral substitution sites to tetrahedral substitution sites was set as 0:6, 2:4, 3:3, 4:2 and 6:0, and was denoted as 0-6, 2-4, 3-3, 4-2 and 6-0, respectively. Eventually, twenty types of Na-MT crystal models with different crystal chemistry properties were established. For example, when the layer charge density of Na-MT was 0.38, the ratio of octahedral sites to tetrahedral sites was 0:6, it was recorded as Na-MT10-6, and more detailed description was shown in Table 1. Natural Na-MT usually contains one layer H2O and the value of c is about 12.5 Å [5,35,36]. To simulate the influence mechanism of crystal chemistry properties on the spatial distribution of H2O and Na + , sixty four H2O molecules were added into the interlayer domain of Na-MT randomly by the adsorption location module to form one layer of H2O, and the the corresponding c value of Na-MT was about 12.5 Å. In the adsorption location module, the parameters were set as follow, the number of cycles was 5, steps per cycle were 5 × 10 4 , geometry optimization and automated temperature control were selected. Three dimensional periodic boundary conditions were applied to all of the simulation systems [17,32]. The established models of Na-MT are shown in Figure 1 (Take Na-MT10-6 and Na-MT16-0 as examples).

Simulation Parameters
The geometry optimization and molecular dynamics simulation were accomplished based on the Forcite module of Material Studio 7.0 [37]. Firstly, the geometry optimization task was performed to fully relax the structure of the Na-MT model under the following conditions. (1) The Na-MT sheets were assumed to be rigid that the lattice parameters of a, b, α, and γ are fixed, but c and β are variable [38]. (2) The algorithm was set to smart, the RMS force was 0.1 kcal/mol/Å, the energy difference was 2 × 10 −5 kcal/mol [39], the RMS displacement was 1 × 10 −5 Å, and the RMS stress was 100 kPa. Secondly, after geometry optimization, the NPT ensembles (isobaric-isothermal ensemble: particle number (N), pressure (P) and temperature (T) are constant.) were run to balance the structure of Na-MT-water systems. The initial velocities were random, temperature was 298.0 K, pressure was 1 × 10 −4 GPa, the time step was 0.5 fs, total simulation time was 500 ps, number of steps were 1 × 10 6 , the temperature control method was Nosé-Hoover [37], and the thermostat was Berendsen. Then, the NVT ensemble (canonical ensemble: particle number (N), volume (V) and temperature (T) are constant, total momentum are also constant) was employed to obtain atomic trajectories of particles in optimized systems [19]. In this part, the other parameters were the same as the last step, but total simulation time was 600 ps, it was run to ensure that the system has reached equilibrium, and the final models were used for analysis [26,40]. During the simulation, the simple point charge (SPC) model was used to represent the interaction of H2O [41,42], in which the geometrical structure of these molecules is rigid, the bond length is 1.0 Å and bond angle is 109.5° [43]. The CLAYFF force field which has been proved to work very well in previous simulation studies about MT was employed to describe the interatomic potential for Na-MT-water systems [44][45][46][47]. The short range van der Waals interaction was calculated based on the atom-based method, and the long-range electrostatic interaction was calculated by the Ewald summation method [48,49]. The hydration models of Na-MT for different layer charge density characteristics after geometry optimization are shown in Table 2   (a) (b) Figure 2. The models of Na-MT after geometry optimization, (a) Na-MT10-6; (b) Na-MT16-0.

Results and Discussion
To comprehensively study the influence of the crystal chemistry properties on the spatial distribution of particles in the Na-MT interlayer domain, the z-axis concentration profiles, self-diffusion coefficient and radial distribution function (RDF) were analyzed, respectively.

Z-Axis Concentration Profiles of Particles in Interlayer of Na-Montmorillonite (Na-MT)
The z-axis concentration profiles can reflect the density characteristic of particles in interlayer of Na-MT [4], therefore, the density variation tendency of H2O and Na + can be predicted and the influence mechanism can be explained by analyzing the simulation results. In this part, it was calculated on the basis of NVT simulation results and the protocols are set as follows: the number of bins was 100, the width of smoothed profile was 1.0 bin.

Z-Axis Concentration Profiles of H2O
The z-axis concentration curves of H2O in interlayer of Na-MT with different crystal chemistry properties are shown in Figures 3 and 4.
It can be seen from Figure 3 that H2O are mainly distributed near the central axis of interlayer domain in the range of −2.0 Å-+2.0 Å, but with the increase of layer charge density from 0.38 to 0.70, the relative concentration of H2O decreases from 10.20 to 8.55, which indicates that H2O move gradually towards both sides of Na-MT sheets with the increase of layer charge density. The reason is that with the increase of layer charge density of Na-MT, the polarity of Na-MT surfaces is strengthened, and the hydrogen bonding force between H2O and tetrahedral surfaces also increases. Therefore, the relative concentration of H2O decreases with the increase of layer charge. This conclusion fits well with our previous study [29].
Compared with the effect of layer charge density, the impact of isomorphic substitution on the distribution of H2O in Na-MT is limited, and it can be seen from Figure 4 that when the layer charge of Na-MT is low, the relative concentration curves of H2O almost coincide, such as Na-MT1 and Na-MT2. But when the layer charge of Na-MT is high, with the amount of octahedral substitution sites increasing, the strength of corresponding peaks decreases gradually, such as Na-MT37-2 and Na-MT39-0. This phenomenon is more obvious in Na-MT4. The reason is that octahedral charges locate in the middle of Na-MT layers, which are farther away from H2O than that of tetrahedral charges. When octahedral substitution sites increase in Na-MT with high layer charge, the hydrogen bonding force between both sides of Na-MT sheets and HW (H in H2O) are weakened because of the long-range spatial effect, which results in H2O no longer being trapped in the central axis of the interlayer space [43]. This simulation result is markedly different from our previous study about Ca-MT that isomorphic substitution has no significant effect on the z-axis density distribution of H2O [39]. It may be attributed to the lower basal spacing of Na-MT than that of Ca-MT, the distance between Ot (O in tetrahedron) and HW is closer in Na-MT than that of Ca-MT, so the impact of isomorphic substitution on the distribution of H2O is more obvious in Na-MT.

Z-Axis Concentration Profiles of Na +
The z-axis concentration curves of Na + in interlayer of Na-MT with different crystal chemistry properties are presented in Figure 5.
As shown in Figure 5a, we can find that layer charge density has an obvious impact on the distribution of Na + in Na-MT. The Na + are mainly distributed near the central axis of the interlayer domain in the range of −1.5 Å-+1.5 Å, but with the increase of charge density from 0.38 to 0.70, the relative concentration of Na + decreases clearly from 12.0 to 7.0. This simulation results show that Na + tend to locate near the central axis of interlayer space in Na-MT1, and diffuse to both sides of Na-MT sheets with the increase of layer charge. The reason is that with the increase of the layer charge density of Na-MT, the negative charge sites generated by isomorphic substitution of Na-MT increase. Furthermore, under the action of electrostatic gravity, Na + are attracted to the vicinity of the tetrahedral surfaces. As a result, the relative concentration of Na + decreases. This result is consistent with our previous study on the distribution of Ca 2+ between Ca-MT layers [39].
It can be seen from Figure 5b that the relative concentration curves of Na + are unimodal and change gradually to irregular distribution as the amount of tetrahedral substitution sites increases. This indicates that Na + are distributed along the central axis when isomorphic substitution mainly occurs in octahedral layers of Na-MT, but Na + tend to move towards tetrahedral sheets when the amount of tetrahedral substitution sites increases [4]. In Na-MT with the same layer charge, Na + are closer to the tetrahedral sheets than octahedral layers, when the amount of tetrahedral isomorphic substitution sites increases, the electrostatic interaction between Na + and Na-MT tetrahedral isomorphic substitution sites will be strengthened due to the short-range spatial effect. Therefore, Na + will shift from the central axis towards both sides of tetrahedral sheets with the increase of tetrahedral substitution sites.
(a) (b) Figure 5. The z-axis concentration curves of Na + in interlayer of Na-MT, (a) effect of layer charge density on the distribution of Na + ; (b) effect of isomorphic substitution sites on the distribution of Na + in interlayer of Na-MT4 (the concentration distribution of Na + in interlayer of Na-MT1, Na-MT2 and Na-MT3 is similar to that of Na-MT4, thus here only Na-MT4 is displayed.)

The Self-Diffusion Coefficient of Particles in Interlayer of Na-MT
The particles in the interlayer of Na-MT will move according to the equations of motion that define the system and have a tendency away from the original position when the system reaches equilibrium [20]. The diffusion coefficient, which can be represented by one-sixth of the slope of mean square displacement (MSD) curves, reflects the change degree of the interlayer particles with time. In order to study the mobility of H2O and Na + , a three-dimensional Einstein relation was used to calculate the self-diffusion coefficient: Here | − 0 | represents the MSD of interlayer species, N is the total number of a given species in the simulation cell, ri(t) is the center-of-mass position of the ith atom at time t [43]. For diffusion coefficient, the results are gotten including two steps, the first step is that the MSD curves are outputted from the last 100 ps of NVT simulation, the second step is that the diffusion coefficient is obtained by linear fitting the MSD curves using Origin 9.0, and the one-sixth of the slope of regression equation is diffusion coefficient. Finally, the results are shown in Table 3 and Figure 6. Table 3. The self-diffusion coefficient of H2O and Na + in interlayer of Na-MT with different crystal chemistry properties (unit: ×10 −10 m 2 /s, and the standard error were generated from linear fitting in Origin 9.0).

Species
Type 0   It can be seen from Table 3 and Figure 6a that the isomorphic substitution of Na-MT can obviously influence the mobility of H2O, with the increase of the octahedral substitution ratio, the self-diffusion coefficient of H2O in the interlayer of the four types of Na-MT increases. The reason is that the attraction force between Na-MT sheets and H2O is weakened because of the long-range spatial effect when the ratio of octahedral substitution increases. Moreover, it can also be seen from Figure 6a that the self-diffusion coefficient of H2O decreases with the increase of layer charge density. The reason is that with the increase of layer charge density, the available free space within the Na-MT interlayer is decreased (Table 2). At the same time, the polarity of Na-MT surfaces is enhanced, and the attraction force between H2O and Na-MT sheets increases, which results in the decrease of the self-diffusion coefficient of H2O [20].
It can also be seen from Table 3 and Figure 6b that the effect of crystal chemistry properties of Na-MT on the self-diffusion coefficient of Na + is similar to that of H2O. With the increase of the tetrahedral substitution ratio or the layer charge density of Na-MT, the self-diffusion coefficient of Na + decreases. When the layer charge density of Na-MT is 0.70 and the ratio of octahedral substitution sites to tetrahedral substitution sites is 0, the self-diffusion coefficient of Na + and H2O reaches the minimum value of 0.0940 × 10 −10 m 2 /s and 0.6381 × 10 −10 m 2 /s, respectively, compared with H2O, the diffusion coefficient of Na + is obviously lower. It can be explained that H2O is neutral and Na + is positive, thus Na + are fixed more easier than that of H2O by negatively charged Na-MT surfaces.

Radial Distribution Function (RDF) of Particles in Interlayer of Na-MT
The RDF is defined as the probability density of finding atom A at a distance away from atom B [46], which plays an important role in describing the aggregation characteristics between two types of atoms. The equation for RDF is: Here is the radial distribution for atom B around atom A, is the number density of atom B, r is the distance between atom A and atom B, NA-B is often denoted as coordination number N(r) for B around A. In this work, the RDF of Ot-Hw and Ot-Na, as well as the Na-Ow were investigated (Hw and Ow stand for the atom H and atom O in H2O respectively, Ot stands for the atom O in tetrahedral surfaces of Na-MT). In this part, the protocols are set as follows: the cutoff radius is 10.0 Å, the interval is 0.02 Å, and periodic self-interactions are included.

The RDF of Ot-Hw
The RDF of Ot-Hw in the interlayer of the four types of Na-MT with different crystal chemistry properties is shown in Figures 7 and 9, respectively. It can be seen from Figure 7 that the first peak of g(Ot-Hw) in Na-MT with different charge density appears at a distance of approximately 1.676 Å. Moreover, the intensity of the first peak of g(Ot-Hw) increases clearly from 0.5 to 0.8 with the increase of layer charge density. The result indicates that H2O are more tightly packed together near the tetrahedral surfaces in the interlayer domain of Na-MT. The reason is that with the increase of the layer charge density of Na-MT, the polarity of Na-MT surfaces increases, therefore, H2O are more likely attracted by stronger hydrogen bonding force and eventually form hydrogen bonds with Ot at the distance of 1.676 ± 0.043 Å away from tetrahedral surfaces (Figure 8, Na-MT40-11 is shown) [29,30].  It can be seen from Figure 9 that when the layer charge density of Na-MT is lower (a, b and c), the isomorphic substitution has no obvious effect on the g(Ot-Hw); only when the layer charge density is higher (d), the g(Ot-Hw) increase with the increase of amount of tetrahedral substitution sites. This phenomenon is consistent with the z-axis concentration distribution of H2O in the Section 3.1.1.

The RDF of Ot-Na
The influence of the crystal chemistry properties on g(Ot-Na) in Na-MT is shown in Figure 10 and Figure 11.
As can be seen from Figure 10, with the increase of layer charge density, the strength of the first peak of g(Ot-Na) increases. The reason is that the higher the layer charge density of the Na-MT, the stronger the polarity of the Na-MT surfaces. The simulation result indicates that under the action of strong polar force, most of the Na + are attracted to the Na-MT surfaces and coordinate with Ot in the tetrahedral surfaces at a position of 2.38 Å. Compared with the coordination distance of 1.676 ± 0.043 Å between Hw and Ot which can be seen in Figure 8 in Section 3.3.1, the coordination distance is farther.
As can be seen from Figure 11, the impact of isomorphic substitution on the strength of g(Ot-Na) is slight when the layer charge density of Na-MT is lower (a, b and c). But when the layer charge density is higher, the strength of the first peak decreases obviously (d, Na-MT411-0). The reason is that tetrahedral charges are closer to Na + than that of octahedral charges in Na-MT layers. The electrostatic attraction between both sides of Na-MT sheets and Na + becomes weakened because of the long-range spatial effect when the octahedral substitution ratio increases, which results in the strength of the first peak of g(Ot-Na) decreasing in Na-MT411-0.

The RDF of Na-Ow
The RDF of Na-Ow in the interlayer of Na-MT with different layer charge density is presented in Figure 12. It can be seen from Figure 12 that the layer charge density of Na-MT has a significant effect on the strength of g(Na-Ow), the main peak of g(Na-Ow) appears at a distance of approximately 2.38 Å and the distance does not change any more with the increase of the layer charge density, which shows that H2O are fixed near the Na + with a Van der Waals radius of 2.38 Å in Na-MT interlayer. The pictorial diagram about the coordination number and Van der Waals radius of Na + are shown in Figure 13. This research result is in line with the conclusions acquired by Sun et al. and Seppälä et al. [6,50]. Besides, the strength of the main peak of the g(Na-Ow) increases gradually with the decrease of layer charge density. The reason is that the higher the charge density, the stronger the polarity of Na-MT surfaces. Most H2O in Na-MT interlayer are adsorbed on the Na-MT surfaces under the action of strong polarity gravity, and form hydrogen bonds with the Ot in the tetrahedral surfaces at a position of 1.676 ± 0.043 Å away from the tetrahedral surfaces ( Figure 8). Therefore, the coordination strength between Na + and atoms O in H2O is weakened. Figure 14 is the RDF of Na-Ow with different isomorphic substitution sites. It can be seen that the isomorphic substitution sites in the structure of Na-MT can significantly influence the aggregation characteristics of Na-Ow ion pairs. With the increase of the octahedral substitution ratio, the strength of the main peak of g(Na-Ow) increases gradually. The reason is that with the increase of the amount of octahedral substitution sites, the distance between the negative potential and Na + increases, and the electrostatic attraction between Ot in the tetrahedral surfaces and Na + in Na-MT interlayer decreases due to the long-range spatial effect. As a result, fewer Na + are transferred to the both sides of tetrahedral surfaces of Na-MT, and more Na + coordinate with Ow in H2O near the central axis, which leads to the enhancement of the strength of g(Na-Ow).

Conclusion
In this work, the influence of layer charge density and isomorphic substitution sites on the spatial distribution of H2O and Na + in the interlayer of Na-MT was investigated. The simulation results show that the distribution of H2O and Na + is significantly related to the crystal chemistry properties of Na-MT. In the interlayer domain of Na-MT, H2O are mainly distributed in the range of −2.0 Å~+2.0 Å near the central axis, and when layer charge density increases, H2O tend to shift and form hydrogen bonds with Ot in tetrahedral surfaces at a distance of 1.676 ± 0.043 Å. The impact of isomorphic substitution on the distribution of H2O depends largely on the charge density of Na-MT, when layer charge density is high (such as Na-MT3 and Na-MT4), H2O in the interlayer domain obviously move towards both sides of Na-MT sheets with the increase of octahedral substitution ratio. Compared with H2O, Na + are mainly concentrated in the range of −1.5 Å~+1.5 Å near the central axis of the Na-MT interlayer, and when layer charge density increases, Na + also move towards Na-MT sheets and coordinate with Ot at a distance of 2.38 Å. But when layer charge is fixed, the diffusion trend of Na + is opposite to that of H2O, and Na + diffuse obviously to both sides of Na-MT sheets with the increase of tetrahedral substitution ratio. The diffusion coefficients of both H2O and Na + decrease when layer charge density or tetrahedral substitution ratio increases. The radial distribution characteristics of Na-Ow reflect that the aggregation degree of Na-Ow ion pairs becomes worse with the increase of layer charge density or tetrahedral substitution ratio, and Na + are hydrated by four H2O molecules at a Van der Waals radius of 2.386 ± 0.004 Å.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, The convergence plot of the system, calculation method and standard error of diffusion coefficient. Funding: This work was supported by grants from the National Nature Science Foundation of China: "Design of structure and gelling performance of Montmorillonite/Alkylammonium based on the adsorption properties of Alkylammonium on Montmorillonite" (No: 51774200) and "Study on the immobilization effect and mechanism of heavy metal in sulfide-mine tailings by chelating agent modified bentonite" (No: 51764003).