Molecular simulation in surface hydration of clay minerals: a review of theory and applications

Clay minerals, which are prevalent gangue minerals, are found in tailings and beneficiation effluent after the extraction of valuable minerals. The surface of clay mineral particles is easy to hydrate, which makes it the main factor restricting tailings separation and wastewater treatment. However, the microscopic mechanism of clay mineral particle surface hydration is not yet systematic. In recent years, with the development of molecular simulation theory and the improvement of computational efficiency, density functional theory (DFT) and molecular dynamics (MD) have gradually become a powerful tool for studying the surface hydration of clay mineral particles, which provides new insight into the interaction between the crystal structures of clay minerals and the interfacial interaction in surface hydration of clay mineral particles at the molecular or atomic levels. This article first reviews the basic theory of DFT and MD, then reviews the research progress on clay mineral surface hydration. From the perspective of molecular simulation, a comprehensive discussion of the clay mineral phase structure, the establishment of the supercell surface model, the clay-water interface interaction and the limitations of molecular simulation was conducted. Water molecules can adsorb with different mineral surfaces in slime water through hydrogen bond, which is the basis of surface hydration mechanism. The hydration layer is composed of three water layers with different densities, with a thickness of about 8-10 Å. Water and ions form hydrate cations, which are adsorbed on the surface of clay minerals, change the structure of water layer on the surface of minerals. This article ends with a brief discussion of conclusions and perspectives. Page 2 of Min et al. Miner Miner Mater 2022;2:3 https://dx.doi.org/10.20517/mmm.2022.01 18


INTRODUCTION
Clay minerals, most of which are layered silicate minerals, can be usually divided into two-layer structure unit layers according to the structure unit layer, namely 1:1 type (such as kaolinite minerals, etc.) and threelayer structure unit layer, namely 2:1 type (such as mica and smectite minerals, etc.) [1] . There are also other different classifications, for example, Ismadji et al. [2] classified clay minerals into four groups: kaolinite, illite, smectite, and vermiculite. Due to the unique interface properties, clay minerals are very important and widely used in various industries [3,4] . Clay minerals with commercial value include kaolinite, montmorillonite, and illite, which are the most prevalent types [5] . Moreover, a large amount of clay minerals is in industrial wastewater, such as mineral separation, papermaking and sewage treatment [6,7] . The clay mineral surfaces have the characteristics of strong hydrophilicity, strong electronegativity, large specific surface area, and are easy to expansion and sludge. Especially in the mineral processing industry, the presence of clay minerals has a significant impact on the separation of useful minerals, the treatment of tailings and beneficiation wastewater, and surface hydration is the dominant factor. Correspondingly, the wide application of clay minerals increases the difficulty and complexity of water treatment to a certain extent.
Currently, the most often utilized techniques for studying clay minerals mainly include surface force apparatus [8] , atomic force microscopy [9,10] , viscosity method [11,12] , and X-ray reflectivity [13,14] . Apart from this, some other methods of indirectly measuring the amount of water adsorbed on the surface of clay mineral particles. According to Einstein's viscosity equation, Min et al. [15,16] introduced the hydration index I as a means of determining the degree of hydration of small clay mineral particles in various electrolyte aqueous solutions. Additionally, they examined the influence of sodium and calcium electrolytes on the hydration degree of montmorillonite and kaolinite aqueous dispersion using the hydration index I. Although there are many related studies on the theoretical analysis and determination methods of clay mineral surface hydration, the understanding of the hydration mechanism is still not clear enough. Therefore, it is particularly important to analyze the hydration mechanism of clay minerals from the molecular/atomic level.
In recent years, scientists have begun to study the molecular/atomic level adsorption process of water molecules, reagent molecules, and other adsorbates on the surface of clay mineral particles using molecular simulation approaches. However, there is no systematic review article on the simulation calculation method of clay mineral surface hydration process and the analysis of hydration mechanism. This article mainly discusses in detail from the basic theory of molecular simulation and the method of building models. Especially, the molecular simulation methods and applications for studying the hydration mechanism of clay mineral particles are systematically organized and reviewed. Finally, we also analyzed the limitations of molecular simulation in model construction, force field selection and other parameter settings.

Density functional theory
Over the past 30 years, density functional theory (DFT) has become a frequently used tool to study the physical and chemical properties of materials. The basis of DFT has been described in many articles and books [17][18][19] . In 1926, Schrödinger defined the rule of particle motion using a linear partial differential equation [20] . Solving Schrödinger equation can be regarded as the basic problem of quantum mechanics.
However, the solution of the Schrödinger equation is a complex multi-body problem. DFT intuitively originated from an approximate electronic structure theory proposed by Thomas [21] , Fermi [22] , and Dirac [23] . They proposed to simulate the kinetic and exchange energy of electronic systems through the uniform electron gas energy density, but the theory is too rough to bind molecules [24] . Hartree-Fock theory uses the Slater determinant to approximate the N-electron wave function, because the Slater determinant satisfies the conditions of the Pauli Exclusion Principle of electrons. Although the Post-Hartree-Fock approach is widely established and has a high degree of accuracy, it is unsuitable for calculating this sophisticated multielectronic system. In 1964, Hohenberg-Kohn theorems [25] establish the feasibility of obtaining ground state properties by varying the electron density, but the exchange-correlation function (E XC [{Ψ i }]) has not yet been defined. In 1965, Kohn and Sham solved this problem [26] . They regarded the interaction of all electrons as an electronic system, and merged this electron with its own Coulomb interaction (Self-Interaction) into the exchange-correlation function for modification. Until this time, the density functional theory was fully established.
Kohn and Sham proposed a simple model for E XC , the so-called "local density approximation". In fact, the charge density is not uniform, and a class of functionals including charge density and density gradient is developed. This is known as the generalized gradient functional (GGA) [27] . And to describe the exchange energy using a mixed form of the exact exchange energy and the GGA exchange functional, such as the Becke 3-parameter-Lee-Yang-Parr (B3LYP) function [28] . Only when the exchange-correlation function used is determined can a DFT be used to describe a specific problem. The most widely used functional calculation of clay minerals is GGA-PBE. The interaction between water and clay minerals is studied at the GGA-PBE level through DFT. From the frontier molecular orbits energy gap, adsorption energy, Mulliken charge transfer, electron density difference analysis and electronic state density to analyze the adsorption mechanism of water and clay minerals [29][30][31][32][33]1] . The B3LYP functional has been proved to be a good predictor of proton transfer processes and changes in free energy in water systems. This functional has also been used in clay mineral surface hydration studies [34][35][36][37] .

Molecular dynamics
Although DFT solves the problem of quantum mechanics, its computational power increases exponentially with the number of electrons, so it is still difficult to describe the quantum mechanical motion and chemical reaction of a larger system. An important calculation method, molecular mechanics (MM) has been developed. MM is based on classical mechanics theory and uses empirical potential functions to characterize the interaction between structural elements; by solving Newton's equations of motion, the trajectories of solid phase points are depicted to screen out energy extreme points and corresponding molecular conformations, calculate the equilibrium or non-equilibrium properties. The MM calculation method does not include electron motion, instead, the system's potential energy is determined by computing the position of the atomic nucleus, which simplifies the calculation.
Molecular simulation approaches include the molecular dynamics (MD) and Monte Carlo (MC) methods, respectively. MD enables the calculation of energy or interactions between molecules by changing the system's molecular coordinates. MC calculates the energy between different conformations the difference completes the change of conformation. These two simulation techniques are widely used in various fields, and the research system decides which simulation method to use. Generally speaking, when the system contains dynamics properties, molecular dynamics is the best choice, MC is more suitable for non-timedependent statistical average calculation. The advantages and disadvantages of molecular simulation calculations are directly controlled by the reliability and efficiency of the potential function.
The interaction potential between molecules described by the mathematical analytical formula of the empirical potential function is called the molecular force field. Defining a force field requires not only determines the function form, but also determines a large number of force field parameters. Various types of force fields have been developed, for example, OPLS [38] , AMBER [39] , CHARMM [40] , CVFF [41] , COMPASS [42] , UFF [43] , DREIGING [44] . These force fields are capable of reproducing a wide variety of properties of carbohydrates, solvents, surfactants, polymers and proteins, DNA, drug molecules are quite consistent with the experiment. Apart from the force fields described above, there are certain customized force fields that have been developed for particular systems on the basis of common force fields. In the research field of clay minerals, the commonly used force fields are ClayFF and polymer consistent force field (PCFF). Cygan et al. [45] developed the ClayFF force field to simulate hydration and multi-component mineral systems, as well as their interactions with aqueous solutions. This method mainly relies on the interaction of van der Waals and Coulomb to maintain the coordination configuration of the central metal cation of the tetrahedron or octahedron. Heinz et al. [46,47] developed PCFF-phyllosilicates and PCFF-Interface force fields for inorganic-organic and inorganic-biomolecular interface based on PCFF. Chen et al. [1] have used ClayFF to study the structure of the water layer on the kaolinite surface, the results show that there are three ordered hydrated films on the kaolinite surface.

Bulk and surface models of clay mineral
The unit cells of clay minerals are established mostly based on unit cell experimental parameters. Researchers have obtained unit cell parameters of clay minerals through techniques such as X-ray and neutron powder diffraction. For example, Bish et al. [48,49] used low temperature (1.5 K) neutron and CuKα Xray powder diffraction to obtain the kaolinite unit cell parameters, respectively. Lage et al. [50] studied the effect of water and alcohol on the desorption of kaolinite from fused pyridinic heterocycles using the unit cell and structure of kaolinite from the article of Bish. Mellini and Zanazzi [51] . provided the crystal structural parameters of lizardite-1T and lizardite-2H1. Demichelis et al. [52] established the unit cell of lizardite and compared it with the calculation model. The montmorillonite model used by Sun et al. [53] in MD simulation also came from the unit cell parameters measured by X-ray powder diffraction [54] . Muscovite unit cell with Si→Al substitution pattern according to Ref. [29] . Si NMR was used to developed supercell (5 × 3 × 1) model [55,56] . Materials Studio provided a unit cell model of clay minerals, and after initial relaxation, the relaxed structure's properties accord with experimental data [57][58][59] .
To produce a layer surface, expand and cleave the unit cell model of clay minerals. Creating a surface of an acceptable size is to accommodate the adsorbate. Constructing a vacuum layer along to the c-axis is to avoid contact between neighboring slabs. The thickness of the vacuum layer is usually 1-2 times that of the surface model in the DFT computation. The thickness of the vacuum layer is 70-80 Å in MD simulation, which may be effective to eliminate the influence of periodic boundary conditions [56] .
DFT has been widely used in simulating electronic properties and reactivity of crystals, surfaces, and molecules. The DFT method is also often used to study the hydration [31,1] , adsorption [60][61][62] , plasticity [63] , electronic [64] and mechanical properties [65] of clay minerals in the last two decades. Clay minerals, composed of aluminum-oxygen octahedral and silicon-oxygen tetrahedral sheet structures. They're divided into "layer types" based on the ratio of tetrahedral and octahedral sheets that have joined together, for example, kaolinite is a typical 1:1-type, the typical 2:1-type dioctahedral clay minerals such as montmorillonite, illite and pyrophyllite [66] . Unit structure of kaolinite based on the experimental value [67] , and Chen [68] has verified the accuracy of kaolinite model. Figure 1A, it has the alumina octahedral surface (AlO 6 octahedral layer) and chemically inert silica tetrahedral surface (SiO 4 tetrahedral layer), corresponding to (001) surface [shown as Figure 1C] and () surface [shown as Figure 1D], respectively [70] . The surface configurations can be cut from the kaolinite supercell shown Figure 1B. The physicochemical properties and utilization of clay minerals are related to their interface structure and properties.

As shown in
In fact, in the process of mineral formation, mining and sorting, due to the formation of environment, thermal vibration, external stress or radiation, the lattice points inside the mineral lattice will have a certain degree of deviation, leading to the existence of lattice defects in the crystal structure of minerals. Therefore, in addition to the ideal mineral model, the adsorption characteristics of minerals with lattice defects can also be studied in clay mineral correlation simulation. Chen et al. [71] analyzed the type and content of Fe in kaolinite from Huaibei mining area by Mossbauer spectrometer. Then, the lattice substitutions of different valence Fe in kaolinite crystal were calculated by DFT, and a reasonable supercell model of kaolinite surface was established. The lattice impurities of Fe 2+ /Fe 3+ in kaolinite and their influence on the surface adsorption performance of kaolinite were studied. Man et al. [72] constructed undoped, Mg, Ca, and Fe (II)-doped kaolinite (001) surfaces, all of which consist of plates of six atomic layers and a 20 Å vacuum zone. All calculations were performed in a four-formula unit, Al 2 Si 2 O 5 (OH) 4 supercell, an aluminum ion is replaced by a defective ion. The effects of Mg, Ca, Fe (II) doping on the atomic structure of kaolinite (001) surface and the subsequent adsorption and infiltration of H 2 O between layers have been systematically studied. Voora et al. [73] is a montmorillonite model based on pyrophyllite model, which is constructed by atomic substitution and adding balanced cations between layers. The balanced cations are generally alkali ions or alkaline earth ions. Therefore, the molecular formula of M-MMT supercell with M as cation is M 2 , where x = 2, M represents a univalent basic cation, x = 1, M represents a bivalent basic cation. Lavikainen et al. [74] designs the model by calculating the position of magnesium replacement in the octahedral layered structure. The addition of the first magnesium substitution was chosen arbitrarily, and the second substitution was selected at different positions. By detecting the total energy of the system, the molecular formula of the final substitution model was determined as Si 32 Al 14 Mg 2 O 80 (OH) 16 . Yl et al. [75] split the surface of smithsonite (101) based on the optimized smithsonite crystals and used it to establish the surface model. At the same time, 15 Å vacuum layers were constructed. Water was placed in 15 × 15 × 15 Å cells for optimization. It allows the water and the top two atomic layers to relax, while the bottom three are fixed. The results show that the adsorption of water on the surface of perfect smithsonite (101) is stronger and more stable.

DFT calculations of H 2 O on clay mineral surfaces
The hydration of clay mineral refers to the state and ability of adsorbing molecules on the surface of clay mineral particles in mud. The surface of clay mineral particles can adsorb polar water molecules directly. More importantly, the surface of clay mineral particles with negative charge, can absorb a large number of cations dispersed in the mud to achieve hydration. It is the structure water in hydration layer. It is the key factor of agglomeration settlement and dehydration in water treatment process. In order to explore the mechanism of mineral surface hydration, so as to effectively control its hydration properties. The adsorption behavior of water molecules on clay mineral surfaces has been studied extensively.
The active sites of the mineral surface should be determined firstly. It can be explored by comparing the adsorption energies of adsorbates at potential adsorption sites of the adsorbent. The optimal active site has the largest adsorption energy, and the optimal adsorption configuration corresponds to it. This method also can be used in searching the active site on the unknown surface. Water molecule adsorption on the active sites of the clay mineral surface has been carried out [31] . Water molecules interact with the montmorillonite basal and edge surfaces are shown as (A) and (B) in the Figure 2, and the interaction of water molecules on Reprinted from Ref. [69] , copyright (2018), with permission from Elsevier. reprinted from Ref. [31] , copyright (2016), with permission from Elsevier. (C) is kaolinite (001) surface; (D) is kaolinite (00 ) surface; the adsorption sites are denoted as T (top positions), B (bridge positions), and H (hollow positions) in the picture; reprinted from Ref. [1] , copyright (2019), with permission from Elsevier.] the Al-OH and Si-O surfaces of kaolinite are shown as (C) and (D) in Figure 2. Comparing and analyzing the adsorption energy of each potential site, and the optimal adsorption configurations of water molecule adsorption on the montmorillonite surfaces are shown as Figure 3.
Adsorption of water molecules on various kinds and surfaces of montmorillonite was investigated. Water molecule adsorption on the Li-Mt surface enhances Li + diffusion [76] . Electron density would be transmitted from Na-Mt to water molecules. Water has a higher adsorption energy on the (010) edge surface than on the (001) basal surface. And water adsorption on the Na-Mt (001) and (010) surfaces were dominated by hydrogen bonding and electrostatic interaction, respectively [31] . Du studied water adsorption on the ammonium illite (001) surface with lattice substitution and the ammonium illite () surface without lattice substitution. The results show that the optimal active site of the illite (001) surface is located at the oxygen atom bonded to the substituting atom Al, with 2 hydrogen bonds forming and the adsorption energy is -0.67eV. As to the () surface, the active region mainly located above the silicon-oxygen ring holes. The adsorption energy is -0.41eV, and water molecules and surface oxygen atoms form three hydrogen bonds. Hydrogen bonds and electrostatic interactions are important in adsorption [77] . Similarity, Zhang investigated the adsorption of water molecules on the Al-OH and Si-O surfaces of kaolinite. It appears that the Al-OH/water system is more stable than the Si-O/water system. According to the findings. The strongest hydrogen bonds are those between water molecules, and those between water molecules and Al-OH are greater than those between Si-O and water [78] . Chen estimated the water molecule's adsorption energies on the potential adsorption sites on the kaolinite (001) and () surfaces. It is found that the adsorption energies of water molecule are -72.12 to -19.23 kJ/mol and -19.23 to -5.77 kJ/mol, respectively. Water molecules preferentially adsorb on the kaolinite (001) surface, with electron transfer between interacting atoms and the Mulliken charge population, showing that the water molecules are adsorbed on the kaolinite surface via hydrogen bonding [1] . The interaction between surface active sites and water molecules, which is usually a combination of hydrogen bonding and electrostatic contact, is the basic hydration process of the clay mineral surface.

MD simulation in surface hydration characterization of clay mineral
Compared with the DFT method, the MD method can simulate the adsorption status of multilayer water molecules or other adsorbents on the mineral surface. For example, clay mineral hydration is bulk water molecules adsorption onto the mineral surface.
At atmospherically relevant temperatures of 298 K and 235 K, structure and water uptake on the kaolinite basal and edge surfaces were investigated. The water models TIP5P-E and SPC/E were taken into consideration. Water molecules easily form monolayers on the Al surface, but not on the Si surface. The water density ρ(z) in the structure reached a maximum of 3-5 g/cm 3 , showing that the monolayers are dense [79] . To realize the efficient clarification and dewatering of high muddy slime water, Min's research group simulated and revealed the hydration properties of kaolinite, montmorillonite, and illite. As the Figure 4 shows that on the kaolinite surface, different monolayer coverage rate (ML) based on the surface oxygen atom (Os) number have the various adsorption configurations. In terms of the quantity of hydrogen bonds produced, water cluster adsorption on the kaolinite (001) surface was more powerful than that on the kaolinite (00_1) surface.
The concentration distribution profiles of Ow (O atom of water) in the Z direction of kaolinite surface show the distribution of the water molecules [shown as Figure 5]. Indicating water molecules adsorption on the kaolinite surfaces, forming a hydration layer composed of three layers of water molecules with a thickness of 8-10 Å. Beside kaolinite, the simulation results of montmorillonite and illite also confirmed the hydration of reprinted from Ref. [31] , copyright (2016), with permission from Elsevier. (C) is kaolinite (001) surface; (D) is kaolinite (00 ) surface; reprinted from Ref. [1] , copyright (2019), with permission from Elsevier.] clay mineral surfaces.
Minerals not only interact with water molecules, but the solution properties also affect the hydration characteristics of the mineral surface. The structure, properties and thickness of the hydration layer are affected by a series of factors. Such as solvated ions and the solution pH. Cation adsorption on mineral surfaces produces a synergistic effect with water molecules. The outer-sphere and inner-sphere are the possible adsorption styles. Generally, outer-sphere complex indicate the weakly interaction, some water molecules exist in the region between clay mineral surface and cation; while the inner-sphere complex indicates that the cation interacts significantly with the surface [80] .
Cs + , Na + , Cd 2+ , and Pb 2+ adsorption on the kaolinite basal surface was investigated using the MD technique at concentrations of 0.1, 0.5, and 1.0 M. Chloride ions formed inner-sphere adsorption at the aluminum vacancy on kaolinite (001) surface, and chloride ions promoted cesium ions and sodium ions to form innersphere adsorption on kaolinite surface at high concentration. Cadmium and lead do not form inner-sphere complexes. On the siloxane surface of ditrigonal cavities, relatively strong outer-sphere cadmium and lead Reprinted from Ref. [1] , copyright (2019), with permission from Elsevier. ML: Monolayer coverage rate.  [1] , copyright (2019), with permission from Elsevier. complexes can be found [81] . Illite is highly selective for adsorption of Cs. When Cs is adsorbed on the illite surface, it will combine with other alkali cations, and the overall selective order is Li + > Na + > K + > Rb + . However, the percentage of Cs assigned to inner-sphere complexes as determined by NMR is in the sequence Na + > K + > Rb + > Li + , which differs from the order of Cs selectivity, and the divergence grows as the molar fraction of Cs in solution decreases [82] . Surface charge properties of kaolinite are affected by divalent ion adsorption, and strong hydrating ions have a significant impact on dispersion pH. With the kaolinite surface, the Fe 2+ ion form an inner spherical complex, whereas the Ca 2+ and Mg 2+ ions form an outer spherical complex [83] . Apart from the basal surface, ions can also be adsorbed by the interlayer surface and pores of clay minerals. The structure and transport of aqueous solutions in 6 nm thick clay mineral nanopores were explored, as well as the effect of clay mineral composition on model Na-hectorite and Namontmorillonite surfaces. And the results demonstrate that the solution structure and diffusion properties are similar at each surface, with minor changes in Na + adsorption complexes and water structure in the first adsorbed layer by reason of various layer hydroxyl group configurations in the two clay mineral models. Surprisingly, surface disturbance affects bulk-like solution structure and diffusion only in a few water layers [84] . Montmorillonite is a clay mineral with a three-layer sheet structure composed of aluminum oxide octahedron in the middle and silicon oxide tetrahedron in the upper and lower layers. Water may be adsorbed in the interlayer and the structural hydroxyl groups' cis-vacancy arrangement in the octahedral is predicted to have a major impact on the adsorbed water's structure. The tetrahedral sheet forms a higher number of hydrogen bonds with bridging oxygens, whereas the octahedral sheet forms weak hydrogen bonds with water in the clay mineral interlayers. Water adsorbed in the interlayer with higher ordering and lower density appears to be more "ice-like". It elucidates the energetics of water adsorption, which determines the swelling, transport, and fracture characteristics of water on a pore-to-macroscopic scale [85] . The wettability of montmorillonite (001) surfaces was evaluated using MD simulations as a function of surface charge and exchangeable cations. Water molecules have a low adsorption energy on the negatively charged montmorillonite surface, while exchangeable cations have a significant electrostatic attraction. The capacity of exchangeable cations to produce hydrated cations with multilayer hydrating shells is quite powerful. The ions adsorb as hydration ions on the montmorillonite surface, increasing the clay minerals' total hydrophilicity. As a result, the layer charge is more important to determine surface wettability, and the adsorbed exchangeable cations on the montmorillonite surface is a secondary factor [86] . The interaction of the hydrated ion [UO 2 (H 2 O) 5 ] 2+ with montmorillonite was predicted, and the results suggest that the uranyl aqua ion formed outer-sphere complexes with the clay layers. The hydrate is tightly bonded, with the initial shell and clay oxygens making 1.4 hydrogen bonds. Uranyl-clay interaction sites were discovered as three-Mg substitution groups. Due to partial surface coverage, increasing the uranyl content improves mobility [87] . The surface hydration of clay minerals is a comprehensive process of multi-component. Surface active atoms, solvated ions, water molecules and pH condition affect the hydration characteristics of mineral surface. Fully understanding the hydration properties of clay mineral surface is the theoretical basis for the regulation and high-value utilization of mineral surface.

LIMITATIONS OF MOLECULAR SIMULATION
Although molecular simulation is a strong tool for study at the atomic and molecular levels, and can be used as a complement to experiment and theory, it does have some limits that should be addressed when determining the proper application of this technique. The following subsections provide a quick overview of the major limitations of molecular simulation technology.

Computational model idealization
In molecular simulations, an accurate and practical model determines the accuracy of the calculation results. Unfortunately, the existing calculations of clay minerals used by scholars are idealized in both DFT and MD simulations. The idealization is embodied in the idealization of clay mineral models and hydration environment.
At present, the reported DFT simulations are mainly based on the surface/interface models constructed by the ideal unit cell of clay minerals, which have no deviation [33,[88][89] . Indeed, there are impurity defects in the clay minerals formed under natural conditions, which affect their physicochemical properties and adsorption ability [90] . Kaolinite and montmorillonite are the most studied clay minerals. Taking the example of kaolinite, the perfect bulk phase does not contain any impurity atoms. In order to get closer to the real kaolinite minerals, researchers have carried out lattice substitution in perfect phases, such as replacing Al in the octahedron with Fe, Mg, or Mn, etc. [91][92][93] . On the other hand, in DFT simulation calculation, the number of adsorbed water molecules is limited, which is far from the real clay-water system. Only the active sites of the clay mineral surfaces can be indirectly obtained through the adsorption of a few water molecules [31] .
Similarly, the clay mineral models used in MD to simulate the performance of clay minerals in aqueous environments is either undoped or only partially doped with interlayer ions or impurity atoms. Liao et al. [94] used pure kaolinite as the skeleton to study the decomposition behavior of hydrate in the pure water flow system by MD method. Subramanian and Lammers [95] investigated the thermodynamics of ion exchange (Na + and K + ) focusing on a water-layer hydration state via MD simulations. The reason why MD simulations cannot calculate the doping of multiple atoms simultaneously is that the force field used in MD is difficult to apply to different ions at the same time.

Force field imperfections in MD simulations
The Born-Oppenheimer approximation confirms that the molecular energy can be approximated as a function of the spatial coordinates of the atoms that make up the molecule, including intramolecular bond interactions, intermolecular van der Waals interactions, and intermolecular hydrogen bond interactions for associative systems. Different types of molecular force field describe different molecular systems. The energy function of the force field determines the molecular configuration, and the parameters of the force field determine the interaction energy between molecular groups. The force field functions are derived mainly from empirical formulae of experimental results, so the force field provides inherent approximations of the true forces acting on and between molecules in the simulated system [96] .
Some of these force fields have high precision and strong specificity but are not universal for a certain system. For example, CVFF [41] and CFF [97] force fields have developed from simulating straight-chain alkanes and cycloalkanes to simulating molecular complexes and crystals. AMBER [39] and CHARMM [40] fields are mainly used to simulate biomolecular systems such as proteins. MM [98] field has excellent performance in reproducing molecular potential energy surface and calculating molecular gas phase properties. The OPLS [38] field, developed by Jorgensen's group at Yale University, focuses on the simulation of organic and polypeptide molecular systems. Some force fields were built to simulate more systems, such as the COMPASS [42] and DREIDING [44] field developed by Sun et al. [42] , which covered almost all elements of the periodic table, but the accuracy of calculation was reduced.
The ClayFF and PCFF fields are mainly used in MD simulations of clay mineral surface hydration. Cygan et al. [45] proposed a universal force field, ClayFF, which has been widely used to simulating hydration and multi-component mineral systems, as well as their interactions with aqueous solutions. The ClayFF force field is protable because it treats the majority of interatomic interactions as non-bonds while maintaining a limited set of parameters that allows for modeling of relatively large and highly disordered systems. However, it should be noted that the ClayFF field focuses on the field parameters of clay minerals, which can support the study of water molecules on the mineral surface. If a more complex collector system containing organic molecules is studied, the ClayFF field does not have the field parameters of organic molecules, so other fields with better description of organic molecules need to replace. Argyris et al. [99] used ClayFF field to simulate the structure and dynamics of water on the surface of graphite, quartz and alumina, and set three kinds of quartz with different degrees of surface basification in the simulation. It is found that solid surface roughness, chemical heterogeneity and surface charge heterogeneity affect the structure and dynamic process of water at the interface, and gives the preferred adsorption sites of oxygen atoms at the interface.
PCFF field belongs to the second generation of force field, based on FF91 field, suitable for polymer and organic materials. The PCFF field is more accurate because of the consideration of many cross-acting items. Heinz et al. [46,47] took nano organic composite materials as the research object and proposed PCFF phyllosilicates' force field with layered silicate potential energy parameters based on PCFF field. Due to the accurate calculation of PCFF for organic molecules, reasonable and consistent simulation results can be obtained for the adsorption behavior of organic molecules at the most important solid-liquid interface in flotation system.
Even though there are many force fields today, no force field can cover all molecules due to the diversity of molecular forms. Therefore, when dealing with practical applications, especially when dealing with multimolecule systems, it is often impossible to find all corresponding parameters in the same force field, that is, the phenomenon of missing parameters. Although not all parameters can be found in one force field, the parameters required by these molecules can often be obtained from multiple other force fields, so a combination of force fields is necessary. For example, Du and Miller [100] used the ClayFF-Dreiding field to describe the interaction of Dodecyl Trimethyl Ammonium Bromide with dextrin molecules on talc surfaces. Simulation results show that the hydrophobic talc base plane does not closely contact with water molecules, and there is a gap of size 3Å on the base plane. Wang et al. [101] combined the ClayFF field with the Dreiding field, which has the ability to predict the domain of organic molecules, to explore the adsorption characteristics and mechanisms of amine surfactants on the surface of fossil quartz solvent with a wide range of pH values by MD simulation. Wang et al. [102] used the Clayff-Dreiding field to study the molecular dynamics of a DDA (DDA) -quartz system under different pH conditions. In recent years, researchers have made an effort to develop ReaxFF force field, which has the potential to dramatically improve the prediction capacities of simulations for a variety of complicated and critical processes [103,104] .

Limitations of simulation time and size scales
As a key parameter affecting the integral process of molecular dynamics, the increase of integration step can accelerate the calculation speed, but too large integration step will cause the divergence of integration process and the inaccuracy of integration value. The main limitation of choosing integral step is the high frequency motion of the system [105,106] . In molecular dynamics simulation, the more calculation steps are, the longer the time step of a single step is, and the closer the final simulation results are to the real situation, which represents the reasonable sampling period of the simulated macro phenomenon. At present, most MD simulations are carried out between picoseconds and microseconds [1,107] . However, due to the limitation of calculation conditions, the simulation time cannot be long enough. In addition, a longer time step will lead to too frequent collisions between particles, system data overflow; smaller time steps reduce the ability to search for phase space during simulation. Therefore, it is necessary to select the appropriate time step according to the actual situation of the simulation system for simulation calculation [108][109][110] .
In molecular dynamics simulation, the scale limitation of the simulation system will cause a certain degree of boundary effect, and periodic boundary conditions can eliminate this effect. Periodic boundary conditions refer to the fact that there are countless periodic replicas around the simulation box, while in the actual calculation process, only N atoms in the main box are explicitly considered. Because this is done under periodic boundary conditions, as soon as one of the atoms leaves the main box, the atom passes back from the other side. After introducing the boundary conditions, some properties of the system will be changed, making the simulation results different from the real system. Therefore, the influence of periodic boundary conditions on the system needs to be corrected within a certain range [111][112][113] . And for different shapes of simulation boxes, there are different correction coefficients when performing diffusion calculations [114] . For convenience, cube simulation box is generally used for simulation calculation [115] . Jamali et al. [116] studied the scale effect of molecular diffusion in cube and cuboid simulation boxes, respectively.
The results showed that the diffusion coefficient perpendicular to the long side of the cuboid simulation box increased with the increase of the side length, and the diffusion coefficient along the long side was basically unchanged or slightly decreased, but this phenomenon was not found in the cube simulation box. Kikugawa [117] found that the diffusion coefficient generated by the simulation is always less than its volume value by establishing a cuboid simulation box. It is found that the periodic boundary conditions have different corrections to these non-cubic simulation boxes in different directions.
Yang et al. [118] confirmed this view to a certain extent. Rozmanov found that the diffusion rate in the corresponding area direction was proportional to the square root of the area when he used the rectangular simulation box with different areas in three directions for calculation. There are also some other polyhedrons that can fill the whole space without gaps that can be used as simulation boxes for simulation calculation. Among them, the most prominent is the truncated octahedral simulation box and the rhombic regular dodecahedron simulation box [119,120] . Because they are more similar to the sphere than the cube or cuboid simulation box, it is of special significance for the simulation study of isotropic fluid structure. At the same time, the distance between the periodic images of these two simulation boxes is larger, so they are more suitable for the study of macromolecules. For example, Papavasileiou et al. [121] studied the interaction between DNA molecules and fullerenes; Ghadari et al. [122] studied the binding characteristics between paclitaxel derivatives and proteins by molecular dynamics simulation.

CONCLUSIONS AND PERSPECTIVES
The clay mineral surface hydration system has been reviewed extensively using DFT and MD simulation calculations, including the interaction between the clay surface and water, the arrangement and structure of the water layer on the clay surface, and the synergistic effect of ions or hydrated ions and water molecules on the clay surface. The development and theory of DFT and MD, as well as the structure of clay minerals, were also briefly discussed. It provides a foundation for larger-scale and high-complexity model simulations as molecular simulation theory and supercomputing efficiency methods advance.
Molecular simulation shows that the adsorption capacity of water and the sodium montmorillonite edge surface is stronger than that of the base surface, and the electron density is transferred from the surface to the water molecule. Water is easily adsorbed on the substituted atoms of illite, hydrogen bonds and electrostatic interaction dominate the entire adsorption process. On the kaolinite surface, water preferentially adsorbs on the (001) surface of kaolinite and is more stable, driven by hydrogen bonding. The adsorption of water and clay mineral surface is the basis of the surface hydration mechanism. Multiple layers of water form a stable hydration structure on the surface of the clay mineral. The hydration layer is composed of three water layers with different densities, with a thickness of 8-10Å. The doping, defect and substitution of clay minerals or exist in an ionic solution environment. The surface of clay minerals interacts with ions, these ions and water molecules are adsorbed on the surface of the minerals in concert. Water and ions form hydrated cations, the adsorption of the inner spherical complex on the clay mineral surface is stronger than that of the outer spherical complex. Hydrated cations are adsorbed on the surface of clay minerals, changing the structure of the water layer on the surface of the minerals, enhancing the hydrophilic and hydrating properties of the clay minerals.
Depending on previous extensive investigations of simulation calculations and experimental studies, the molecular simulation of clay mineral surface hydration has the following conclusions and application prospects: (1) The mechanism of adsorption of different dissociated surface water molecules in slime water was comprehensively and systematically studied to master the surface hydration microscopic mechanism of other mineral particles and real slime particles in slime water except clay mineral particles.
(2) Aiming at the control mechanism of interface hydration of coal slime particles, research was carried out to optimize the molecular structure design of hydrophobic modified agents by weakening hydration film, so as to provide theoretical and technical support for the development of new agents.
(3) Clay mineral surface hydration is a complex process influenced by a variety of circumstances, there is no perfect force field for surface hydration. Through the combination of high-precision hydration test and simulation, the force field used to describe the hydration of clay mineral surface is comprehensively developed, particularly incorporating the interaction of hydrated ion and clay mineral surface or interface.

Authors' contributions
Made substantial contributions to conception and design of the study and supervised, as well as validated: Min F Collated and analyzed data, wrote original draft, reviewed and edited: Wang L Conducted investigation and supervision: Chen J, Liu C Performed data acquisition and validation: Ren B, Zhang L, Zhu Y

Availability of data and materials
Not applicable.