Melting Behaviour under Pressure of Kaolinite Clay: A Nanoscale Study

: In this study, the curves of variation of melting temperature as a function of pressure were determined for pressures up to 20 GPa using molecular dynamics (MD) calculations. The CLAYFF force ﬁeld is used for the simulated PT curve of the clay kaolinite structure. For this purpose, we have adopted the Z-method to determine the melting point ( Tm ) and superheat limit temperature (TLS) for different densities in kaolinite clay. In addition, various quantities, such as the radial distribution function (RDF), the mean square displacement (MSD), and the diffusion coefﬁcient were evaluated in order to ensure the solid behaviour at the superheat limit temperature and the liquid behaviour at the melting point for the equilibrated structure of kaolinite.


Introduction
Nanoclays consist of one or two-dimensional octahedral and tetrahedral sheets stacked together and can be classified into two types based on their arrangement of T and O layers: tetrahedral-octahedral (TO) and tetrahedral-octahedral-tetrahedral (TOT).Examples of TO nanoclays include kaolinite, halloysite, and serpentine, while montmorillonite, hectorite, synthetic laponite, sepiolite, and palygorskite fall under TOT fully or partially stacked structures [1][2][3].Kaolinite is an uncharged layer of clay with a dioctahedral phyllosilicate structure.It has a 1:1 layered aluminosilicate.The layers stack along the c-axis and are composed of a repeating layer of an octahedral sheet of aluminium (O) and a tetrahedral sheet of silicon (T) with a basal distance from the elemental layer which varies from 7.1 to 7.4 Å.Each layer is formed of a sheet of tetrahedral SiO 4 forming six-membered silicate rings connected by common oxygen atoms to a sheet of octahedral AlO 6 forming fourmembered aluminate rings.The structure of kaolinite clay and its transformations during heating have been studied for many years [4][5][6][7][8][9][10][11].However, a huge lack of information regarding the knowledge of thermodynamic properties with or without external stress still exists.
The thermodynamic properties of kaolinite under high-temperature conditions have not been sufficiently studied.Understanding the pressure dependence of melting temperature is crucial for a solid material and especially for research fields related to geology and earth sciences.Nowadays, molecular dynamics (MD) simulation is a well-established technique in various fields, including soil and mineral sciences.This is particularly true for the case of studies related to high pressures and high temperatures, which can hardly be achieved experimentally.Among the superheating and melting studies based on MD calculations are those of Matsui and Price, Chaplot et al., and Belonoshko et al. [12][13][14][15].
Murray presented the melting point of kaolinite in pure state at 2123 K under atmospheric pressure [5], but in a clayey form, whereas Benazzouz et al. [9] used the molecular Minerals 2023, 13, 679 2 of 14 dynamics simulation calculation to evaluate the melting point at 1818 K at higher pressure conditions equal to 8.85 GPa.This melting point value was investigated using the Z-method and was used to investigate the phase diagram of kaolinite [9].From computer simulations of melting, it is observed that the limit superheating (T LS ) is higher than the melting temperature (T m ) by about 20%-30% [15].
In this work, the structural bulk of kaolinite under a high-temperature environment was used to evaluate the melting temperature (T m ) and the limit of superheating (T LS ).In order to do so, we used the developed Z-method, which was widely applied in materials modelling.The theory of the Z-method was first established by Belonoshko et al. [15].The Z-method helps to evaluate more easily the superheating and melting temperature values, which are difficult to achieve experimentally.Some applications using the Z-method were already devoted to various studies such as the properties of the Lennard-Jones FCC crystal model at the limit of overheating [16], molybdenum at high pressure and temperature [17], the atomistic model for homogeneous fusion [18], high-pressure-hightemperature polymorphism in Ta [19], melting temperature of kaolinite [9], melting curve of silica at high pressures [20], and the phase diagram and equation of state of chromium [21].
The MD simulations described here use the CLAYFF force field developed by Cygan et al. to simulate the kaolinite structure [22].One of the main features of the CLAYFF force field is its flexibility within the clay lattice where the metal-oxygen interactions are described by a van der Walls term and a Coulomb function.First, the CLAYFF force field has already proven to be very effective in modelling the crystal structure of different phases of hydroxide, oxyhydroxide, and clay [22].Then, it was used for different materials such as the structural, dynamic, and energetic properties of cementitious materials [23], the study of the structural and vibrational properties of talc and pyrophyllite [24], the thermomechanical properties of the Montmorillonite crystal, insulated clay nanoplate [25], the simulation of hydration, and the elastic properties of montmorillonite [26].Moreover, this method was used in other previous topics to investigate, for instance, layered and nanoporous materials and their aqueous interfaces [27].Based on the above-mentioned studies, we investigate here the high-pressure melting of kaolinite using the molecular dynamics method.
The paper is organised as follows.Section 2 gives the computational details and method related to this study.Section 3 presents the results and discussion.The conclusion is given in the final section.

Model and Computational Details
The atomic structure of kaolinite clay is the same as that described in our previous study [9].The structure of kaolinite shown in Figure 1 is originally based on crystallographic data from Bish et al. [28], a unit cell of triclinic symmetry with lattice parameters a = 5.153 Å, b = 8.941 Å, c = 7.390 Å, and angles: α = 91.926, which corresponds to 34 atoms, and its density is 2.6 g/cm 3 .
The simulation box used in our calculations is 4 × 2 × 3 unit cells in dimensions a, b, and c (24 kaolinite unit cells), with a total of 816 atoms in the solid.The corresponding Lx, Ly, and Lz dimensions are 20.614Å, 17.881 Å, and 22.739 Å. Periodic boundary conditions are applied in all three directions.
The potential energy of a simulated system includes the non-bonded energy and the bonded energy.The non-bonded energy components are Coulomb interactions and short-range van der Waals energy.The potential energy of a simulated system includes the non-bonded energy and the bonded energy.The non-bonded energy components are Coulomb interactions and shortrange van der Waals energy.
The first non-bonded energy term (Coulomb interactions) varies just like the inverse of the bond length between the i and j ions (rij).The following equation shows this energy: where qi and qj represent the charge of ion i and j respectively, derived from quantum mechanics calculations, e is the electron charge, and εo is the dielectric permittivity of vacuum (8.85419 × 10 −12 F/m).
The second non-bonded energy term (van der Waals interactions) presents the interaction between the two non-bonded atoms, represented by the conventional Lennard-Jones (12-6) function, and it includes the short-range repulsion associated with the increase in energy as the two atoms approach each other (Urep), the attractive term energy, (Uattr), and the London dispersion forces: where Do and Ro are the empirical parameters derived from the fitting of the model; this is equivalent to the dissociation energy and equilibrium atomic separation, respectively.The parameters of van der Waals between the interaction of two atoms are calculated in accordance with the arithmetic and the geometric mean rules for the distance parameter, Ro, and the energy parameter, Do, respectively [29].The first non-bonded energy term (Coulomb interactions) varies just like the inverse of the bond length between the i and j ions (r ij ).The following equation shows this energy: where q i and q j represent the charge of ion i and j respectively, derived from quantum mechanics calculations, e is the electron charge, and ε o is the dielectric permittivity of vacuum (8.85419 × 10 −12 F/m).
The second non-bonded energy term (van der Waals interactions) presents the interaction between the two non-bonded atoms, represented by the conventional Lennard-Jones (12-6) function, and it includes the short-range repulsion associated with the increase in energy as the two atoms approach each other (U rep ), the attractive term energy, (U attr ), and the London dispersion forces: where Do and Ro are the empirical parameters derived from the fitting of the model; this is equivalent to the dissociation energy and equilibrium atomic separation, respectively.The parameters of van der Waals between the interaction of two atoms are calculated in accordance with the arithmetic and the geometric mean rules for the distance parameter, Ro, and the energy parameter, Do, respectively [29].
Minerals 2023, 13, 679 The values of the non-bonded interaction parameters are given in Table 1.The bonded energy is only used to describe hydroxyl groups.The bond stretch energy of the hydroxyl bond (O-H) can be described by a harmonic (simple quadratic) expression: where k 1 is an empirical force constant, r o is the equilibrium bond length, r is the atomatom separation distance.The values for the CLAYFF-bonded parameters for hydroxyl interactions are given in Table 2. Molecular dynamics calculations at the thermodynamic equilibrium were performed using Nosé-Hoover method in the DL_POLY code, developed in Daresbury Laboratory [30].
The ensembles used in most MD simulations are isothermal-isobaric (NPT), microcanonical (NVE), canonical (NVT), and grand-canonical (µVT), where the variables N, P, T, V, E, and µ represent, respectively, the number of particles, the pressure, the temperature, the volume, the energy, and the chemical potential.In this work, the standard ensembles implemented are the NPT (the number of particles N, the pressure P, and the temperature T are fixed) with barostat and thermostat relaxation times of 0.5 and 0.5 ps, respectively; and, in the second calculation, the NVE (the number of particles N, the volume V, and the energy E are fixed during the simulation).The Nosé-Hoover thermostat was used to control the temperature [31].The pressure was kept constant using a Melchionna modification of the Hoover algorithm [32] in which the equations of motion involve a Nosé-Hoover thermostat and a barostat.A thermostat and a barostat were used for the NPT ensemble.
The motion equations were integrated by using the Verlet leapfrog integration algorithm [33,34] with a time step of 1 fs (1 × 10 −15 s) and 0.5 fs to ensure the conservation of energy.The Ewald summation technique [35] is considered as the most satisfactory method for treating these Coulombic long-range interactions, using a tolerance of 1 × 10 −5 .For short-range interactions, either spherical or "minimum image" cutoff criteria are commonly used [29,34].The cutoff radius r c = L x /2, where L x is the smallest box length (L x = 1.788 nm).For van der Waals interactions, the cutoff used is between 0.8 and 0.85 nm.
Different samples (solid structure) were prepared at 200 picosecond (ps) using equilibration NPT MD simulation and performed at atmospheric pressure 0.1 MPa [9], 1 GPa, 5 GPa, and 10 GPa at 300 K. Figure 2 shows the kaolinite structure at 0.1 Mpa after equilibration.

Z-Method to Quantify Melting
The Z-method [15,16] is used to find the melting temperature ™ and superheating limit temperature (TLS) of kaolinite clay.The temperature of melting occurs when the Gibbs free energies of the solid and liquid phases are equal, as shown in the calculations of Morris et al. [36].Unlike the two-phase (also called co-existence) method, molecular dynamics simulations, the Z-method uses single-phase (solid).The Z-method takes its name from the characteristic shape of the isochoric curve.In this work, the calculations have been achieved in the NVE (N number of atoms, volume V, energy E) ensemble.The volume remains fixed throughout the isochoric curve.
The isochore consists of three branches (See the figure below with the Z curve).The first is crystalline (the solid branch).It can be seen that, before arriving at the temperature limit of superheating (TLS), the structure remained solid because the system of interaction of atoms is equilibrating when the initial kinetic energy is low, and then as this energy starts to increase with T, the system approaches and enters in the superheated solid phase across the melting curve.With all these different temperatures, we obtain an almost straight line to the point of TLS, which is the upper cap of the letter Z.
A second branch with negative slope is a transition from a crystalline state to a liquid state.Once TLS isochore is reached, the temperature drops to Tm due to the latent heat of melting (i.e., the energy required to change the phase of a substance from solid to liquid) being removed from the kinetic energy.This part is called the intermediate branch that starts with TLS and ends with Tm, which is a discontinuity in the straight line.Finally, once the temperatures are above Tm, the system begins to be liquid and then forms an almost straight line (a liquid branch) that is the lower cap of the letter Z.
Therefore, the point at the end of the solid branch represents the temperature limit of superheating, TLS, and the first point which is at the beginning of the curve is the melting temperature, Tm.All details regarding this method are provided in the works of Belonoshko et al. [15,16].

Results and Discussion
In a previous work, we have determined the melting point and the temperature limit of superheating at pressure, using the Z-method [9].Here, we will determine the melting and superheating limit curves.

Z-Method to Quantify Melting
The Z-method [15,16] is used to find the melting temperature ™ and superheating limit temperature (T LS ) of kaolinite clay.The temperature of melting occurs when the Gibbs free energies of the solid and liquid phases are equal, as shown in the calculations of Morris et al. [36].Unlike the two-phase (also called co-existence) method, molecular dynamics simulations, the Z-method uses single-phase (solid).The Z-method takes its name from the characteristic shape of the isochoric curve.In this work, the calculations have been achieved in the NVE (N number of atoms, volume V, energy E) ensemble.The volume remains fixed throughout the isochoric curve.
The isochore consists of three branches (See the figure below with the Z curve).The first is crystalline (the solid branch).It can be seen that, before arriving at the temperature limit of superheating (T LS ), the structure remained solid because the system of interaction of atoms is equilibrating when the initial kinetic energy is low, and then as this energy starts to increase with T, the system approaches and enters in the superheated solid phase across the melting curve.With all these different temperatures, we obtain an almost straight line to the point of T LS , which is the upper cap of the letter Z.
A second branch with negative slope is a transition from a crystalline state to a liquid state.Once T LS isochore is reached, the temperature drops to T m due to the latent heat of melting (i.e., the energy required to change the phase of a substance from solid to liquid) being removed from the kinetic energy.This part is called the intermediate branch that starts with T LS and ends with T m , which is a discontinuity in the straight line.Finally, once the temperatures are above Tm, the system begins to be liquid and then forms an almost straight line (a liquid branch) that is the lower cap of the letter Z.
Therefore, the point at the end of the solid branch represents the temperature limit of superheating, T LS , and the first point which is at the beginning of the curve is the melting temperature, T m .All details regarding this method are provided in the works of Belonoshko et al. [15,16].

Results and Discussion
In a previous work, we have determined the melting point and the temperature limit of superheating at pressure, using the Z-method [9].Here, we will determine the melting and superheating limit curves.In the first time, the sample was equilibrated with the NPT ensemble.We performed simulations to obtain the variation of volume and density of the kaolinite structure as a function of pressure, P, for a temperature T = 300 K.For a given pressure, the points obtained corresponding to volume and density are shown in Figure 3a,b, respectively.We may see that the density of kaolinite structure increases with increasing pressure.
Thereafter, we performed the Z-method simulations, and the computed isochore for each density is shown in Figure 4.This figure shows the melting curve with pressure corresponding to different densities.Each (pressure, temperature) point in the isochoric curve was obtained in the NVE ensemble after performing 800,000 time steps, equivalent to 400 ps.
Minerals 2023, 13, x FOR PEER REVIEW 6 of 14 In the first time, the sample was equilibrated with the NPT ensemble.We performed simulations to obtain the variation of volume and density of the kaolinite structure as a function of pressure, P, for a temperature T = 300 K.For a given pressure, the points obtained corresponding to volume and density are shown in Figure 3a,b, respectively.We may see that the density of kaolinite structure increases with increasing pressure.
Thereafter, we performed the Z-method simulations, and the computed isochore for each density is shown in Figure 4.This figure shows the melting curve with pressure corresponding to different densities.Each (pressure, temperature) point in the isochoric curve was obtained in the NVE ensemble after performing 800,000 time steps, equivalent to 400 ps.In the first time, the sample was equilibrated with the NPT ensemble.We performed simulations to obtain the variation of volume and density of the kaolinite structure as a function of pressure, P, for a temperature T = 300 K.For a given pressure, the points obtained corresponding to volume and density are shown in Figure 3a,b, respectively.We may see that the density of kaolinite structure increases with increasing pressure.
Thereafter, we performed the Z-method simulations, and the computed isochore for each density is shown in Figure 4.This figure shows the melting curve with pressure corresponding to different densities.Each (pressure, temperature) point in the isochoric curve was obtained in the NVE ensemble after performing 800,000 time steps, equivalent to 400 ps.The melting point is estimated at the lowest point of the liquid branch, while the superheated limit temperature is estimated at the highest point in the solid branch.The averaged values of superheated limit temperature and the melting point are shown in Table 3.The obtained melting point values are compared with the unique experimental value given by Murray equal to 2123 K [5].The limit of superheating values is clearly in good agreement with values of temperature superheating calculated using MD calculations in the NVT ensemble given by Benazzouz and Zaoui [8].The superheating temperature simulated are given between 1990 K and 2550 K in pressures ranging from 6 GPa to 20 GPa [8].The limit superheated temperature, T LS , is higher than the melting temperature, T m , in good agreement with a previous work of [15].In Table 3, we present a percentage of this difference.
In order to confirm the behaviour of the solid superheating limit temperature, T LS , and the behaviour of liquid at the melting temperature, Tm, various properties were measured such as the radial distribution function (g (r ij )), coordination numbers (n (r ij )), and mean square displacement (MSD) at different temperatures.The radial distribution function for the various atom-atom pairs is a measure and a good indicator to determine the correlation between atoms within a system.To calculate the RDF between the atoms, we use the following expression: where ρ j is the number density of particles of type j in the system and dn ij is the number of particles of type j situated in a spherical layer of thickness dr at distance r from a particle of type i.The coordination number of particles of type j around a particle of type i gives the value of the following integral.
where ρ j is the number density of the atoms of type j.
Radial distribution functions (RDF) and coordination number, n(r), were calculated for all atom pairs.As an example of this, the Si-O_S (silicon-surface oxygen) and O_H-O_H (oxygen-oxygen) RDF with temperature and their related coordination numbers are presented in Figure 5.In this figure, we show the RDF for density ρ = 2.57 g/m 3 (P = 13.32GPa) at different temperatures 1597.1 K, T LS = 2274 K, 2097 K, T m = 1818 K, and 2390.1 K corresponding to temperatures on the solid branch, the superheating limit temperature, the melting point and the temperature on the branch liquid, respectively.These results are compared to those obtained from the DM calculations using the Z-method [9].
K corresponding to temperatures on the solid branch, the superheating limit temperature, the melting point and the temperature on the branch liquid, respectively.These results are compared to those obtained from the DM calculations using the Z-method [9].
The first peak, corresponding to the shortest distance between two neighbouring atoms, decreases with increasing temperature.This phenomenon is repeated and is clearly perceptible in the figures for the pairs of atoms.It explains the differences in structure and the liquid state of the structure of kaolinite at the melting point.The second properties calculation used to confirm the liquid structure of kaolinite at the melting temperature is mean square displacement (MSD).Figure 6 shows the mean square displacement plotted as function of time (up to 400 ps) for density ρ = 2.76 g/m 3 at 20.07 GPa.This figure presents the MSD of two examples Al and O_H atoms at T = 1726.7K (solid branch), TLS = 2549.5K, Tm = 2290.9K, and T = 2934.7K (liquid branch).The MSD slop at Tm is more important than TLS.It can be noticed that for Al and O_H atoms at T ≥ Tm (liquid branch), the MSD increases rapidly with time, which explains certainly the liquid state of the structure.In the solid branch, the MSD at T = 1726.7K and TLS = 2549.5K increases but with a low value (slopes less weak) for O_H atom.
The diffusion coefficients are calculated to complete the characterisation of the solid and liquid states obtained.The diffusion coefficients of all atoms, Al, Si, O_A (apical oxygen), O_S (surface oxygen), O_H (hydrogen oxygen), and H were computed.The diffusion coefficient, D, is related with the MSD by Einstein expression according to: The first peak, corresponding to the shortest distance between two neighbouring atoms, decreases with increasing temperature.This phenomenon is repeated and is clearly perceptible in the figures for the pairs of atoms.It explains the differences in structure and the liquid state of the structure of kaolinite at the melting point.
The second properties calculation used to confirm the liquid structure of kaolinite at the melting temperature is mean square displacement (MSD).Figure 6 shows the mean square displacement plotted as function of time (up to 400 ps) for density ρ = 2.76 g/m 3 at 20.07 GPa.This figure presents the MSD of two examples Al and O_H atoms at T = 1726.7K (solid branch), T LS = 2549.5K, T m = 2290.9K, and T = 2934.7K (liquid branch).The MSD slop at T m is more important than T LS .It can be noticed that for Al and O_H atoms at T ≥ T m (liquid branch), the MSD increases rapidly with time, which explains certainly the liquid state of the structure.In the solid branch, the MSD at T = 1726.7K and T LS = 2549.5K increases but with a low value (slopes less weak) for O_H atom.
As we may see, the diffusion coefficient is significantly higher in the liquid branch than in the solid branch.A significant increase in the diffusion coefficient is noticed in the liquid isochore branch.Table 4 shows the diffusion coefficient values for all atoms, Al, Si, O_A, O_S O_H, and H at different temperatures for ρ = 2.76 g/m 3 .The value of its diffusion coefficient for Al in the liquid isochore branch obtained from the slope of MSD in Figure 6, D = 4.73 × 10 −9 m 2 /s, is in good agreement with those obtained for the liquid diffusion.The diffusion coefficients are calculated to complete the characterisation of the solid and liquid states obtained.The diffusion coefficients of all atoms, Al, Si, O_A (apical oxygen), O_S (surface oxygen), O_H (hydrogen oxygen), and H were computed.The diffusion coefficient, D, is related with the MSD by Einstein expression according to: where r i (0) and r i (t) are the initial and final positions of the centre of mass of the particle at time t, and <|r i (t) − r i (0)| 2 > is the mean square displacement (MSD) averaged over the ensemble.The diffusion coefficient is calculated graphically from the slope of MSD versus time.
As we may see, the diffusion coefficient is significantly higher in the liquid branch than in the solid branch.A significant increase in the diffusion coefficient is noticed in the liquid isochore branch.Table 4 shows the diffusion coefficient values for all atoms, Al, Si, O_A, O_S O_H, and H at different temperatures for ρ = 2.76 g/m 3 .The value of its diffusion coefficient for Al in the liquid isochore branch obtained from the slope of MSD in Figure 6, D = 4.73 × 10 −9 m 2 /s, is in good agreement with those obtained for the liquid diffusion.Here, it can be seen that the difference in diffusion is clear; the diffusion coefficient is much higher in the liquid branch than in the solid branch.In the solid blanch, the diffusion coefficients are almost zero because the atoms in a solid cannot diffuse away from their equilibrium positions, compared to those of the liquid branch.
Minerals 2023, 13, x FOR PEER REVIEW 10 of 14 Figure 7 presents an example of the variation of diffusion coefficient for Si and O_S atoms as function of temperature for ρ = 2.76 g/m 3 .Here, it can be seen that the difference in diffusion is clear; the diffusion coefficient is much higher in the liquid branch than in the solid branch.In the solid blanch, the diffusion coefficients are almost zero because the atoms in a solid cannot diffuse away from their equilibrium positions, compared to those of the liquid branch.The variation of the diffusion coefficient for O_A and H atoms as function of pressure is shown in Figure 8.The diffusion coefficient decreases with pressure.There is a difference between the values of diffusion at the superheating limit temperature and at the melting point.The diffusion coefficient of these atoms at Tm is greater than the one at TLS by a factor larger than 8.This can be used to confirm once more the liquid phase at Tm.The variation of the diffusion coefficient for O_A and H atoms as function of pressure is shown in Figure 8.The diffusion coefficient decreases with pressure.There is a difference between the values of diffusion at the superheating limit temperature and at the melting point.The diffusion coefficient of these atoms at Tm is greater than the one at T LS by a factor larger than 8.This can be used to confirm once more the liquid phase at T m .
Finally, the points of the melting curve obtained in this work were fitted using a unified equation, proposed by Kechin for the melting curve at high pressure [36,37].The melting equation is given as follows: where P is the pressure, T m (P) is the melting temperature at this pressure, and T 0 is the melting point at zero pressure.The values of the fitted parameters (T 0 , a, b and c) corresponding to melting are given in Table 5. Figure 9 shows the melting curve with the pressure up to 20 GPa and the Kechin fit curve of these points.There is no experimental nor theoretical melting curves to compare with our present work.We may determine the melting points, Tm, at any pressure, the melting point at zero pressure, T 0 , equal to 1100.9 K.The same equation is used to fit the superheating limit temperature curve as shown in Figure 10.The fitted parameters for the Kechin equation are given in Table 5.The variation of the diffusion coefficient for O_A and H atoms as function of pressure is shown in Figure 8.The diffusion coefficient decreases with pressure.There is a difference between the values of diffusion at the superheating limit temperature and at the melting point.The diffusion coefficient of these atoms at Tm is greater than the one at TLS by a factor larger than 8.This can be used to confirm once more the liquid phase at Tm. Finally, the points of the melting curve obtained in this work were fitted using a unified equation, proposed by Kechin for the melting curve at high pressure [36,37].The melting equation is given as follows: where P is the pressure, Tm(P) is the melting temperature at this pressure, and T0 is the melting point at zero pressure.The values of the fitted parameters (T0, a, b and c) corresponding to melting are given in Table 5. Figure 9 shows the melting curve with the pressure up to 20 GPa and the Kechin fit curve of these points.There is no experimental nor theoretical melting curves to compare with our present work.We may determine the melting points, Tm, at any pressure, the melting point at zero pressure, T0, equal to 1100.9 K.The same equation is used to fit the superheating limit temperature curve as shown in Figure 10.The fitted parameters for the Kechin equation are given in Table 5.

Conclusions
In summary, we have calculated the melting point and the superheating limit temperature of kaolinite clay for different densities using the Z-method based on classical molecular dynamics simulation by means of the CLAYFF force field.These important values have been deduced from the simulated isochore (PT) in the microcanonical, NVE, ensemble at different pressures.The obtained melting points, Tm, are 1893 K, 2096.9K, and 2290.9K corresponding to the pressures 9.98 GPa, 13.32 GPa, and 20.07 GPa, respectively.The obtained superheating limit temperatures TLS are 2029 K, 2273.9K, and 2549.5 K related to pressures 7.85 GPa, 12.5 GPa, and 19.15 GPa, respectively.Several measurable quantities, such as the radial distribution function, coordination numbers, mean square displacement, and diffusion coefficients, were used to confirm melting and superheating temperatures in the isochore curve.Finally, the melting and superheating limit curves were plotted for kaolinite clay and their corresponding melting equations were proposed.These melting and superheating limit equations are given as follows:   1100.9 1 /143.6 . .and   1567.4 1 /256.4 . .

Conclusions
In summary, we have calculated the melting point and the superheating limit temperature of kaolinite clay for different densities using the Z-method based on classical molecular dynamics simulation by means of the CLAYFF force field.These important values have been deduced from the simulated isochore (PT) in the microcanonical, NVE, ensemble at different pressures.The obtained melting points, T m , are 1893 K, 2096.9K, and 2290.9K corresponding to the pressures 9.98 GPa, 13.32 GPa, and 20.07 GPa, respectively.The obtained superheating limit temperatures T LS are 2029 K, 2273.9K, and 2549.5 K related to pressures 7.85 GPa, 12.5 GPa, and 19.15 GPa, respectively.Several measurable quantities, such as the radial distribution function, coordination numbers, mean square displacement, and diffusion coefficients, were used to confirm melting and superheating temperatures in the isochore curve.Finally, the melting and superheating limit curves were plotted for kaolinite clay and their corresponding melting equations were proposed.These melting and superheating limit equations are given as follows: T m (P) = 1100.9× 1 + P/143.6) 83.7 e −0.509P and T LS (P) = 1567.4× 1 + P/256.4) 93.4 e −0.326P , respectively.

Figure 1 .
Figure 1.Structure of kaolinite formed by 4 × 2 × 3 unit cells.O atoms are in red, the tetrahedral Si atoms in brown, octahedrally coordinated Al atoms are in grey, and H atoms are white.Projection view along the ZY-direction.

Figure 1 .
Figure 1.Structure of kaolinite formed by 4 × 2 × 3 unit cells.O atoms are in red, the tetrahedral Si atoms in brown, octahedrally coordinated Al atoms are in grey, and H atoms are white.Projection view along the ZY-direction.

Figure 3 .
Figure 3. Variation of the volume (a) and the density (b) as function pressure.

Figure 3 .
Figure 3. Variation of the volume (a) and the density (b) as function pressure.

Figure 3 .
Figure 3. Variation of the volume (a) and the density (b) as function pressure.

Figure 7
Figure 7  presents an example of the variation of diffusion coefficient for Si and O_S atoms as function of temperature for ρ = 2.76 g/m 3 .Here, it can be seen that the difference in diffusion is clear; the diffusion coefficient is much higher in the liquid branch than in the solid branch.In the solid blanch, the diffusion coefficients are almost zero because the atoms in a solid cannot diffuse away from their equilibrium positions, compared to those of the liquid branch.

Figure 7 .
Figure 7. Variation of diffusion coefficient for atoms with increasing of temperature for ρ = 2.76 g/m 3 .

Figure 7 .
Figure 7. Variation of diffusion coefficient for atoms with increasing of temperature for ρ = 2.76 g/m 3 .

Figure 7 .
Figure 7. Variation of diffusion coefficient for atoms with increasing of temperature for ρ = 2.76 g/m 3 .

Figure 8 .
Figure 8. Variation of the diffusion coefficient of O_A and H atoms as function of pressure.

Figure 8 .
Figure 8. Variation of the diffusion coefficient of O_A and H atoms as function of pressure.

Figure 9 .
Figure 9. Melting and Kechin fit curves for kaolinite up to 20 Gpa.Figure 9. Melting and Kechin fit curves for kaolinite up to 20 Gpa.

Figure 9 .
Figure 9. Melting and Kechin fit curves for kaolinite up to 20 Gpa.Figure 9. Melting and Kechin fit curves for kaolinite up to 20 Gpa.

Figure 10 .
Figure 10.Superheating limit temperature and Kechin fit curves for kaolinite with pressure.

Figure 10 .
Figure 10.Superheating limit temperature and Kechin fit curves for kaolinite with pressure.

Table 1 .
Non-bonded parameters for the CLAYFF Force Field.

Table 2 .
Bonded parameters for the CLAYFF Force Field.

Table 3 .
Melting point and superheating limit temperature values for kaolinite with pressure and the percentage between T m and T LS . [9]].

Table 4 .
Values of diffusion coefficient for all atoms, Al, Si, O_A, O_S O_H, and H at different temperatures.

Table 5 .
Fitted parameters for the Kechin equation.

Table 5 .
Fitted parameters for the Kechin equation.