Effect of Layer Charge Characteristics on the Distribution Characteristics of H2O and Ca2+ in Ca-Montmorillonites Interlayer Space: Molecular Dynamics Simulation

The charge characteristics of montmorillonite have significant effects on its hydration and application performances. In this study, a molecular dynamics simulation method was used to study the influence of the charge position and charge density of montmorillonite on the distribution of H2O and Ca2+ in layers. The results showed that when the layer charge is mainly derived from the substitution among ions in the tetrahedron, a large number of Hw and Ot are combined into a hydrogen bond in the interlayer, thus the water molecules are more compactly arranged and the diffusion of water molecules among the layers is reduced. In addition, the Ca2+ are diffused to the sides by a concentrated distribution in the central axis of the layer. As the charge density of the montmorillonite increases, the polarity of the Si–O surface increases, which lesds to the deterioration of the diffusibility of the water molecules and the structure of the water molecules in the interlayers is more stable. The increase in the layer charge density lesds to the expansion of the isomorphic substitution range of the crystal structure, which results in a more dispersed distribution of Ca2+ among the layers under the action of electrostatic attraction between the substituted negative sites and the Ca2+.


Introduction
Montmorillonite is a typical clay mineral with a unique crystal structure, which provides it with a high cation exchange capacity and high expansion capacity, thus, it is widely used in various fields. The molecular simulation method used in this study is useful and has several advantages [1][2][3]. Not only can it simulate the molecular structure of the material [4], it can also simulate the dynamic change of molecules [5,6], intuitively describe the mechanism of the chemical reaction at the molecular and atomic scale, and verify the rationality of conclusions or predict the results of experiments [7,8]. Therefore, many scholars have also used the molecular simulation method to study the hydration properties of montmorillonite [9].
The research study on montmorillonite using molecular simulation technology began in the 1990s. Skipper et al. first used Monte Carlo (MC) and molecular dynamics (MD) methods to simulate the interlayer water of Mg-based montmorillonite (Mg-MMT), and obtained the interlayer spacing of Mg-MMT under different water content conditions [10,11]. When the simulation reached equilibrium, the influence mechanism of montmorillonite layer charge density and inter-ion substitution position on the distribution characteristics of water molecules and Ca 2+ in the interlayer.

Simulation Model Establishment
In this study, Materials Studio 7.0 simulation software was used to simulate the distribution characteristics of Ca 2+ and water molecules in Ca-montmorillonite with different charge density [30]. Material Studio is a material simulation computing platform developed by Accelrys, USA. Firstly, according to the semi-crystal structure formula M x+y (Al 2-x Mg x )(Si 4-y Al y )O 10 (OH) 2 nH 2 O, the corresponding montmorillonite crystal model was established [31]. Where M is the interlayer exchangeable cation (Ca 2+ ) and x + y is the unit layer charge. The structure belongs to the monoclinic C2/m space group, the crystal layer constant a = 0.523 nm, b = 0.906 nm, and the c value is variable, when the structural unit layer is anhydrous, c = 0.960 nm [32]; if water molecules are present among the structural unit layers, the value of c will vary with the number of water molecules and the exchangeable cations in the interlayers. Table 1 shows the atomic coordinates of montmorillonite. Secondly, based on these previously mentioned parameters, the montmorillonite model can be established. The montmorillonite super-cell consisted of 32 formula units in the form of 8 × 4 × 1 in the X-,Y-and Z-axis directions, respectively. Thirdly, according to our studying purpose, we designed different isomorphic substitution conditions in the super-cell. The principle of substitution is: the substitution position is random and two adjacent positions can't be replaced at the same time.
According to the crystal structure characteristics of natural montmorillonite, combined with the research purpose of this paper, the layer charge density of the Ca-montmorillonite semi-cell is 0.25, 0.375, 0.5, 0.625, respectively, and the following were recorded as M1, M2, M3, and M4. The substitution positions occurred in octahedron and tetrahedron, respectively. The ratio of the number of octahedral substitution sites to tetrahedral substitution sites is: 0:6, 2:4, 3:3, 4:2, 6:0, the following are recorded as 0-6, 2-4, 3-3, 4-2, 6-0, eventually twenty montmorillonite crystal models were established. For example, for a montmorillonite with a layer charge density of 0.25, the ratio of the octahedral to tetrahedral substitution sites is 0:6, which is recorded as M1 0-6 .
The water molecules among the TOT (Tetrahedron-Octahedron-Tetrahedron) crystal layers are added by the Adsorption Location module. Since the amount of water molecules in the interlayers can affect the interlayer spacing of the montmorillonite, when the montmorillonite contains a layer of water molecules, the layer spacing is about 12.5 Å, and the two layers are about 15.0 Å [33,34]; in order to study the influence of layer charge characteristics of montmorillonite on water molecules and Ca 2+ in the interlayer, this paper used the established 8a × 4b × 1c super-cell to simulate. Two-hundred-and-fifty-six water molecules were added into the montmorillonite model to form two layers of water molecules, and the corresponding Ca-montmorillonite layer spacing was 15.0 Å [35]. Corresponding to four types of Ca-montmorillonites with different charge densities, five different ion substitution modes and substitution ratios, twenty montmorillonite hydration models were finally established (taking M1 as an example, as shown in Figure 1).

Simulation Parameters
In the geometric optimization, it was assumed that the TOT sheets of the calcium-based montmorillonite remained rigid, the unit cell parameters a, b, α, and γ remained unchanged, and c and β were variables [36,37]. In this study, the Smart Minimizer algorithm was selected. The parameters were set to the following: the RMS Force was 0.1 kcal/mol·Å, the Entergy Difference was 2 × 10 −5 kcal/mol, the RMS Displacement was 1 × 10 −5 Å, and the RMS Stress was 100 kPa. The Ewald summation method was for long-range electrostatic interaction and the atom-based summation method as for the short-range van der Waals forces. The vacuum cutoff spacing was 12.5 Å, the spline width was 1.0 Å, the buffer width was 0.5 Å, and the number of iteration steps was 5000.
The molecular dynamics simulation used the Forcite module and the NPT and NVT ensembles. The temperature was 298.0 K, the step length was 0.5 fs, the total simulation time of the NPT was 100 ps, and the total simulation time of the NVT was 150 ps. The first 50 ps was used for simulation and the last 100 ps was used for analysis. The number of steps was 1 × 10 5 and the output was every 200 steps. The universal force field (UFF) was used for simulation in this work [38]. The UFF is the force field that basically covers all the elements in the periodic table of elements, which can meet the technical requirements for the simulation in this paper. During the simulation, all atoms and ions in the interlayer were unconstrained, allowing atomic coordinates and lattice parameters to change freely, and the output data were used for result analysis. The basic parameters of the twenty Camontmorillonite models are shown in the following Table 2.

Simulation Parameters
In the geometric optimization, it was assumed that the TOT sheets of the calcium-based montmorillonite remained rigid, the unit cell parameters a, b, α, and γ remained unchanged, and c and β were variables [36,37]. In this study, the Smart Minimizer algorithm was selected. The parameters were set to the following: the RMS Force was 0.1 kcal/mol·Å, the Entergy Difference was 2 × 10 −5 kcal/mol, the RMS Displacement was 1 × 10 −5 Å, and the RMS Stress was 100 kPa. The Ewald summation method was for long-range electrostatic interaction and the atom-based summation method as for the short-range van der Waals forces. The vacuum cutoff spacing was 12.5 Å, the spline width was 1.0 Å, the buffer width was 0.5 Å, and the number of iteration steps was 5000.
The molecular dynamics simulation used the Forcite module and the NPT and NVT ensembles. The temperature was 298.0 K, the step length was 0.5 fs, the total simulation time of the NPT was 100 ps, and the total simulation time of the NVT was 150 ps. The first 50 ps was used for simulation and the last 100 ps was used for analysis. The number of steps was 1 × 10 5 and the output was every 200 steps. The universal force field (UFF) was used for simulation in this work [38]. The UFF is the force field that basically covers all the elements in the periodic table of elements, which can meet the technical requirements for the simulation in this paper. During the simulation, all atoms and ions in the interlayer were unconstrained, allowing atomic coordinates and lattice parameters to change freely, and the output data were used for result analysis. The basic parameters of the twenty Ca-montmorillonite models are shown in the following Table 2.

Analysis and Discussion
The distribution characteristics of water molecules and calcium ions in the montmorillonite layer during molecular dynamics simulation can be visually and quantitatively described by the Z-axis density distribution curves, mean square displacement function (MSD), and radial distribution function (RDF) of the final configuration.

Z-Axis Density Distribution
By molecular dynamics simulations, the Z-axis density distribution curves reflects the density distribution of the components in the montmorillonite interlayer in the direction of the vertical montmorillonite layer. Thus, the distribution and relative position of water molecules and Ca 2+ in the interlayer of montmorillonite can be further studied, and the influence mechanism of charge density of montmorillonite on interlayer particles can be explained.

Z-Axis Density Distribution Characteristics of Water Molecules
The Z-axis density distribution curves of water molecules in montmorillonite layers with different layer charge characteristics are shown in Figure 2. It can be seen from the left figure that the Z-axis density distribution of water molecules in the interlayer of Ca-montmorillonite has little difference, and the overall distribution was bimodal. The two peaks were distributed in the range of ±2.5 Å on both sides of the central axis among the two TOT sheets and were substantially symmetrical. Therefore, the different substitution sites had no significant effect on the Z-axis density distribution of water molecules. Figure 2b shows that the layer charge substitution sites had the same characteristics, as the charge density increased from 0.25 to 0.625, The water molecules in the montmorillonite layer move along the axis of the center. Moreover, the water molecules' distribution range extended to both sides along the central axis, the left side moved from −1.7 Å to −2.5 Å, and the right side moved from 1.8 Å to 2.5 Å. The reason was that as the charge density of the montmorillonite layer increased, the polarity of the siloxane surface was strengthened, and the water molecules moved toward the siloxane surface on both sides of the interlayer by the hydrogen bonding force between the oxygen in the siloxane surface and the hydrogen in the water molecules. right side moved from 1.8 Å to 2.5 Å. The reason was that as the charge density of the montmorillonite layer increased, the polarity of the siloxane surface was strengthened, and the water molecules moved toward the siloxane surface on both sides of the interlayer by the hydrogen bonding force between the oxygen in the siloxane surface and the hydrogen in the water molecules.

Z-Axis Density Distribution Characteristics of Ca 2+
The Z-axis density distribution curves of Ca 2+ in the montmorillonite layers with different layer charge characteristics are shown in Figure 3. From the figure we can find that the Ca 2+ in the interlayer were mainly concentrated in the range of −1.5 Å~+1.5 Å along the central axis of the montmorillonite. The change in the charge position in montmorillonite had a certain influence on the distribution of Ca 2+ in the interlayer. When the layer charge as derived from aluminum-substituted silicon in the octahedron, such as Mi 6-0 , the Ca 2+ in the layers mainly showed a single peak and were concentrated in the range of −1.5 Å~+1.5 Å in the interlayer domain. In addition, the different charge densities of the montmorillonite had little difference in the Z-axis density distribution of Ca 2+ , indicating that the layer charge negative potential in the octahedron had little effect on the distribution characteristics of interlayer Ca 2+ because of the long-range spatial effect. As the amount of Mg 2+ substituted Al 3+ in the tetrahedron increases, the amount of Al 3+ substituted Si 4+ in the corresponding octahedron decreases, the Ca 2+ will shift toward the siloxane surface on both sides in the interlayers, and the Z-axis density distribution of Ca 2+ will show a bimodal or trimodal distribution. The phenomenon was more pronounced in the high-charged montmorillonite such as M2, M3, and M4. Due to the short-term effect of space, the Ca 2+ was closer to the negative charge position in the tetrahedron in the interlayer, as part of the Ca 2+ were subjected to a large electrostatic attraction and moved toward the tetrahedron, so the Z-axis distribution of Ca 2+ exhibited a multi-peak distribution.
In addition, the charge density of montmorillonite had a significant effect on the distribution of Ca 2+ in the interlayer from Figure 3. With the increase of the layer charge density, the Ca 2+ moved more significantly from the middle of the layer to the siloxane surface on both sides, as shown in the figure, migration peaks increased or multiple migration peaks appeared. This phenomenon was more pronounced as the proportion of the charge in the tetrahedron increased. The Ca 2+ among the layers were subjected to the negative potential electrostatic attraction of the isomorphism of the tetrahedron, causing the Ca 2+ to diffuse to the sides. Since the montmorillonite with higher charge density had a large number of negative potential points, the Ca 2+ in the interlayer were subjected to strong electrostatic force, and the diffusion to the siloxane surface was more obvious. Therefore, in the same crystal structure of montmorillonite, with the increase of the charge density of montmorillonite, the Ca 2+ appeared to gradually shift from the vicinity of the central axis to the two sides in the montmorillonite layer. large number of negative potential points, the Ca 2+ in the interlayer were subjected to strong electrostatic force, and the diffusion to the siloxane surface was more obvious. Therefore, in the same crystal structure of montmorillonite, with the increase of the charge density of montmorillonite, the Ca 2+ appeared to gradually shift from the vicinity of the central axis to the two sides in the montmorillonite layer.

Mean Square Displacement (MSD) and Self-Diffusion Coefficient
The mean square displacement is the mean square of the change in particle position in respect to its initial position at different times, which can reflect the change of particle offset with time [39]. The diffusion coefficient reflects the mobility of the particles, indicating the magnitude of the change in particle position with time, and its magnitude can be derived from the mean square displacement [40]. By analyzing the diffusion of water molecules and Ca 2+ between Ca-montmorillonite layers, the microscopic mechanism of macroscopic hydration characteristics exhibited by montmorillonite can be studied, which provides theoretical guidance for the deep processing and application of montmorillonite. The self-diffusion coefficient of interlayer particles in the montmorillonite interlayer can be derived from the mean square displacement curve of the particle, and the magnitude of the self-diffusion coefficient is proportional to the slope of the mean square displacement curve, i.e., [41]: where | ( ) − (0)| is the mean square displacement of the molecule, Nα is the total number of α particles, ri denotes the position vector of ith particle and the angular brackets denote an ensemble average, and the self-diffusion coefficient D is 1/6 of the slope of the mean square displacement curve.

MSD and Self-Diffusion Coefficient of Water Molecules
The mean square displacement curves and self-diffusion coefficients of water molecules in the interlayer of the Ca-montmorillonite with different charge characteristics are shown in Table 3 and Figure 4. We can find that the charge position and density of montmorillonite have a certain influence on the diffusibility of water molecules in the interlayer. With the proportion of charge density produced by montmorillonite octahedral substitution increased, the mobility of water molecules in the inter-layer domain was higher and, thus, had a higher diffusion coefficient (Figure 4a, M1). It can be inferred that the self-diffusion coefficient of water molecules was affected by the polarity of the montmorillonite siloxane surface. In other words, if the total charge was fixed, when the charge

Mean Square Displacement (MSD) and Self-Diffusion Coefficient
The mean square displacement is the mean square of the change in particle position in respect to its initial position at different times, which can reflect the change of particle offset with time [39]. The diffusion coefficient reflects the mobility of the particles, indicating the magnitude of the change in particle position with time, and its magnitude can be derived from the mean square displacement [40]. By analyzing the diffusion of water molecules and Ca 2+ between Ca-montmorillonite layers, the microscopic mechanism of macroscopic hydration characteristics exhibited by montmorillonite can be studied, which provides theoretical guidance for the deep processing and application of montmorillonite. The self-diffusion coefficient of interlayer particles in the montmorillonite interlayer can be derived from the mean square displacement curve of the particle, and the magnitude of the self-diffusion coefficient is proportional to the slope of the mean square displacement curve, i.e., [41]: where r i (t) − r j (0) 2 is the mean square displacement of the molecule, N α is the total number of α particles, r i denotes the position vector of ith particle and the angular brackets denote an ensemble average, and the self-diffusion coefficient D is 1/6 of the slope of the mean square displacement curve.

MSD and Self-Diffusion Coefficient of Water Molecules
The mean square displacement curves and self-diffusion coefficients of water molecules in the interlayer of the Ca-montmorillonite with different charge characteristics are shown in Table 3 and Figure 4. We can find that the charge position and density of montmorillonite have a certain influence on the diffusibility of water molecules in the interlayer. With the proportion of charge density produced by montmorillonite octahedral substitution increased, the mobility of water molecules in the inter-layer domain was higher and, thus, had a higher diffusion coefficient (Figure 4a, M1). It can be inferred that the self-diffusion coefficient of water molecules was affected by the polarity of the montmorillonite siloxane surface. In other words, if the total charge was fixed, when the charge source was generated in the tetrahedron, due to the short-range spatial effect, the water molecules were subjected to strong polar attraction and difficult to move, so the self-diffusion coefficient was low.
Compared with the influence of the charge position, the charge density of the montmorillonite had a more significant effect on the diffusibility of the water molecules in the interlayers, and the diffusibility of the water molecules decreased with the increase in the layer charge density (Figure 4b). The reason was that as the charge density of montmorillonite increased, the polar attraction of siloxane to water molecules increased; water molecules moved toward the siloxane surface, and the number of hydrogen bonds formed by the combination of O in the siloxane surface (O t ) and H in the water molecule (H w ) increased. This is consistent with Greathouse's conclusion that water molecular diffusivity is related to hydrogen bond fracture among layers of montmorillonite [3]. Due to the action of hydrogen bonding, water molecules are bound to the surface of the montmorillonite tetrahedron and difficult to diffuse. This phenomenon was more pronounced when the ratio of tetrahedral isomorphism was higher in montmorillonite. number of hydrogen bonds formed by the combination of O in the siloxane surface (Ot) and H in the water molecule (Hw) increased. This is consistent with Greathouse's conclusion that water molecular diffusivity is related to hydrogen bond fracture among layers of montmorillonite [3]. Due to the action of hydrogen bonding, water molecules are bound to the surface of the montmorillonite tetrahedron and difficult to diffuse. This phenomenon was more pronounced when the ratio of tetrahedral isomorphism was higher in montmorillonite.

MSD and Self-Diffusion Coefficient of Ca 2+
The mean square displacement curves and self-diffusion coefficients of Ca 2+ in the interlayer of the Ca-montmorillonite with different charge characteristics are shown in Table 4 and Figure 5. We can see that the self-diffusion coefficients of Ca 2+ were obviously affected by the position of the lattice substitution (Figure 5a). When the layer charge of montmorillonite was mostly formed in a tetrahedron, the diffusion coefficient of Ca 2+ was lower than that formed in an octahedron because of the short-range effect, as mentioned in Section 3.1.2. The influence of charge density on the diffusion of Ca 2+ is shown in Figure 5b; as the charge density of montmorillonite increased, the electrostatic

MSD and Self-Diffusion Coefficient of Ca 2+
The mean square displacement curves and self-diffusion coefficients of Ca 2+ in the interlayer of the Ca-montmorillonite with different charge characteristics are shown in Table 4 and Figure 5. We can see that the self-diffusion coefficients of Ca 2+ were obviously affected by the position of the lattice substitution (Figure 5a). When the layer charge of montmorillonite was mostly formed in a tetrahedron, the diffusion coefficient of Ca 2+ was lower than that formed in an octahedron because of the short-range effect, as mentioned in Section 3.1.2. The influence of charge density on the diffusion of Ca 2+ is shown in Figure 5b; as the charge density of montmorillonite increased, the electrostatic attraction between the Ca 2+ and the Si-O surface strengthened, as a result, the diffusion coefficient of Ca 2+ was lower. This is consistent with Seppälä's [16] conclusion that the diffusion coefficient of calcium ions decreases with the increase in the layer charge.  In summary, Tables 3 and 4 show that the diffusivity of water molecules among montmorillonite layers was much larger than that of Ca 2+ , but the diffusion coefficient of Ca 2+ was more affected by the charge position (Figures 4 and 5), which was because the water molecules were neutral, although the polar side (the side where H was located) would be affected by the negative side of the siloxane surface. However, compared to the Ca 2+ with positive charge, the polarity of the water molecules was much weaker. Thus, water molecules in the interlayer diffused more easily than Ca 2+ . As for the influence from an overall point of view, when the charge density of montmorillonite was high or the ratio of tetrahedral substitution position was high, the hydration performance of montmorillonite was worse due to the low diffusibility of water molecules and Ca 2+ .

The Radial Distribution Function (RDF)
The radial distribution function reflects the probability density among two types of atoms and can reflect the aggregation characteristics of ions in the system [42]. By analyzing the radial distribution function among two different particles in the montmorillonite interlayer domain, we can infer whether there is a force among the two types of atoms and the magnitude of their force. The radial distribution function of the system is the algebraic averaging of the radial distribution functions of all the same atoms, reflecting the atomic distribution characteristics of the system [39]. The calculation formula is: In summary, Tables 3 and 4 show that the diffusivity of water molecules among montmorillonite layers was much larger than that of Ca 2+ , but the diffusion coefficient of Ca 2+ was more affected by the charge position (Figures 4 and 5), which was because the water molecules were neutral, although the polar side (the side where H was located) would be affected by the negative side of the siloxane surface. However, compared to the Ca 2+ with positive charge, the polarity of the water molecules was much weaker. Thus, water molecules in the interlayer diffused more easily than Ca 2+ . As for the influence from an overall point of view, when the charge density of montmorillonite was high or the ratio of tetrahedral substitution position was high, the hydration performance of montmorillonite was worse due to the low diffusibility of water molecules and Ca 2+ .

The Radial Distribution Function (RDF)
The radial distribution function reflects the probability density among two types of atoms and can reflect the aggregation characteristics of ions in the system [42]. By analyzing the radial distribution function among two different particles in the montmorillonite interlayer domain, we can infer whether there is a force among the two types of atoms and the magnitude of their force. The radial distribution function of the system is the algebraic averaging of the radial distribution functions of all the same atoms, reflecting the atomic distribution characteristics of the system [39]. The calculation formula is: where g αβ (r) is the radial distribution of β; n β is the number of β with a radius of r → r + dr; ρ β is the number density of β; and r is the distance between α and β. In this work, the characteristics of the radial distribution function of oxygen in the siloxane surface (O t ) and hydrogen atoms in the water molecules (H w ) and Ca 2+ in the interlayer domain were studied.  Figure 6a shows that first peak of g(O t -H w ) of montmorillonite appears at the distance of approximately 2.5 Å. With the increase of the ratio of isomorphism in the tetrahedron, the intensity of the main peak of g(O t -H w ) was enhanced, which indicates that when the montmorillonite layer charge was mostly derived from the substitution among ions in the tetrahedron, the Si-O surface exhibited a stronger polar force, and the water molecule had a polar side (i.e., the side where H was located), which was subjected to a stronger polar force from the Si-O surface, making the water molecules more compact near the siloxane surface, so the intensity of the main peak of g(O t -H w ) in the RDF was enhanced.
In Figure 6b, the g(O t -H w ) shows the RDF between O t and H w in the interlayer with different charge densities, which indicates that as the layer charge density increased, the intensity of the first main peak of g(O t -H w ) gradually increased. The reason was that the polarity of the Si-O surface increases with the increase in the charge density of the montmorillonite, and the H w was pulled by a stronger polar force when the layer charge was higher; therefore, the first main peak of g(O t -H w ) gradually increased with an increase in the layer charge density. This further demonstrates that the diffusivity of water molecules among the layers of montmorillonite decreased as the layer charge density increased. stronger polar force from the Si-O surface, making the water molecules more compact near the siloxane surface, so the intensity of the main peak of g(Ot-Hw) in the RDF was enhanced. In Figure 6b, the g(Ot-Hw) shows the RDF between Ot and Hw in the interlayer with different charge densities, which indicates that as the layer charge density increased, the intensity of the first main peak of g(Ot-Hw) gradually increased. The reason was that the polarity of the Si-O surface increases with the increase in the charge density of the montmorillonite, and the Hw was pulled by a stronger polar force when the layer charge was higher; therefore, the first main peak of g(Ot-Hw) gradually increased with an increase in the layer charge density. This further demonstrates that the diffusivity of water molecules among the layers of montmorillonite decreased as the layer charge density increased. As can be seen from the left figure, if the layer charge density is fixed, the difference in the substitution positions will significantly affect the distribution of Ca 2+ . When ion substitution occurred, mostly in tetrahedrons (0-6, 2-4, 3-3), Ca 2+ were shifted closer to the Z-axis due to their closer proximity to the substitution source, and the distribution was more uniform. A low peak appeared at a distance (2.8 Å) closer to the Ot, but most of the Ca 2+ were still concentrated in the middle of the layer, about 5.5 Å from the Ot. When the substitution source was mostly in octahedrons (4-2, 6-0), the radial distribution of Ca 2+ was less affected due to the spatial long-range effect As can be seen from the left figure, if the layer charge density is fixed, the difference in the substitution positions will significantly affect the distribution of Ca 2+ . When ion substitution occurred, mostly in tetrahedrons (0-6, 2-4, 3-3), Ca 2+ were shifted closer to the Z-axis due to their closer proximity to the substitution source, and the distribution was more uniform. A low peak appeared at a distance (2.8 Å) closer to the O t , but most of the Ca 2+ were still concentrated in the middle of the layer, about 5.5 Å from the O t . When the substitution source was mostly in octahedrons (4-2, 6-0), the radial distribution of Ca 2+ was less affected due to the spatial long-range effect mentioned in Section 3.1.2.
From Figure 7b, we can find that with the increase in the charge density of montmorillonite, due to the negative charge of the Si-O surface increases, some of the Ca 2+ moved toward the Si-O surface under the action of electrostatic attraction. Therefore, some new peaks appearred around 2.5 Å-3.5 Å, which was more obvious when the layer charge density of the montmorillonite was higher (M3, M4).

Conclusions
By simulating the distribution characteristics of water molecules and Ca 2+ in the interlayer under different charge characteristics of montmorillonite, the following conclusions were obtained: (1) When the layer charge of the Ca-montmorillonite was mainly derived from the replacement of Si 4+ with Al 3+ in the tetrahedron, the Hw and the Ot were more easily combined to form a hydrogen bond in the interlayer. Due to the hydrogen bonding force, the water molecules were more compactly arranged, and the diffusion of water molecules among the layers was reduced. Due to the action of electrostatic attraction between the Ca 2+ and the Si-O surface, the Ca 2+ move to the Si-O surface on both sides and the distribution range was expanded. (2) With the increase of the charge density in the Ca-montmorillonite layer, the polarity of the Si-O surface was enhanced, and the polar force generated by the water molecules in the interlayers was stronger, the number of hydrogen bonds between the Hw and the Ot increased, and the self-diffusion coefficient of the water molecules decreased. For the Ca 2+ , the higher the layer charge density of montmorillonite, the larger the electrostatic force of the Ca 2+ subjected to the negative point, and the larger the distribution range of Ca 2+ .
On the basis of this work, we have a preliminary understanding of the influence of the chargecharacteristic of Ca-montmorillonite on water molecules and Ca 2+ . In the future, we will conduct an in-depth study on the influence mechanism of the charge density of Ca-montmorillonite on the polarity of the Si-O surface [43]. In addition, we will further study the structure, performance, and application characteristics of organo-modified montmorillonite by means of MD. The conclusions about the influence mechanism of layer charge characteristics on the distribution characteristics of H2O and Ca 2+ in Ca-montmorillonites interlayer space will give us useful information about the interlayer of montmorillonite at the atomic level.

Conclusions
By simulating the distribution characteristics of water molecules and Ca 2+ in the interlayer under different charge characteristics of montmorillonite, the following conclusions were obtained: (1) When the layer charge of the Ca-montmorillonite was mainly derived from the replacement of Si 4+ with Al 3+ in the tetrahedron, the H w and the O t were more easily combined to form a hydrogen bond in the interlayer. Due to the hydrogen bonding force, the water molecules were more compactly arranged, and the diffusion of water molecules among the layers was reduced. Due to the action of electrostatic attraction between the Ca 2+ and the Si-O surface, the Ca 2+ move to the Si-O surface on both sides and the distribution range was expanded. (2) With the increase of the charge density in the Ca-montmorillonite layer, the polarity of the Si-O surface was enhanced, and the polar force generated by the water molecules in the interlayers was stronger, the number of hydrogen bonds between the H w and the O t increased, and the self-diffusion coefficient of the water molecules decreased. For the Ca 2+ , the higher the layer charge density of montmorillonite, the larger the electrostatic force of the Ca 2+ subjected to the negative point, and the larger the distribution range of Ca 2+ .
On the basis of this work, we have a preliminary understanding of the influence of the charge-characteristic of Ca-montmorillonite on water molecules and Ca 2+ . In the future, we will conduct an in-depth study on the influence mechanism of the charge density of Ca-montmorillonite on the polarity of the Si-O surface [43]. In addition, we will further study the structure, performance, and application characteristics of organo-modified montmorillonite by means of MD. The conclusions about the influence mechanism of layer charge characteristics on the distribution characteristics of H 2 O and Ca 2+ in Ca-montmorillonites interlayer space will give us useful information about the interlayer of montmorillonite at the atomic level.