# Machine Learning Force Field beyond the Limits of Classical and First-Principles Molecular Dynamics Simulations: The Case of Kaolinite Hydration

David Dell'angelo, Juliette Lainé, Halima Said, Yann Foucaud, Michael

Badawi

<figure><img src="image_1.png" /></figure>

## To cite this version:

David Dell'angelo, Juliette Lainé, Halima Said, Yann Foucaud, Michael Badawi. Machine Learning Force Field beyond the Limits of Classical and First-Principles Molecular Dynamics Simulations: The Case of Kaolinite Hydration. Journal of Physical Chemistry C, 2024, 128 (27), pp.11447-11455. 10.1021/acs.jpcc.4c03288 . hal-04884874

HAL Id: hal-04884874

https://hal.univ-lorraine.fr/hal-04884874v1

Submitted on 13 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Machine learning force field beyond the limits
of classical and first-principles molecular
dynamics simulations: the case of kaolinite
hydration

David Dell'Angelo,*,† Juliette Lainé,‡ Halima Said,† Yann Foucaud,*,¶ and Michael

Badawi*,†

†Université de Lorraine, CNRS, L2CM, F-57000 Metz, France

‡ArcelorMittal Global R&D, 57280 Maizières-lès-Metz, France

¶Université de Lorraine, CNRS, GeoRessources, F-54000 Nancy, France

E-mail: david.dell-angelo@univ-lorraine.fr; yann.foucaud@univ-lorraine.fr;

michael.badawi@univ-lorraine.fr

Abstract

The understanding of the interaction between mineral surfaces and water holds great
significance in unraveling the physical chemistry of surfaces in beneficiation processes.
Within this study, we employ theoretical simulations to explore the adsorption mech-
anism of water molecules on kaolinite surface, a clay mineral commonly found in iron
ores. Albeit ab initio simulations offer key insights into the microscopic aspects of water
adsorption, such as dissociation mechanisms or the hydrophilic/hydrophobic nature of
the surface, the analysis of water structuring is limited by the insufficient statistical
sampling beyond the interfacial monolayer. By extending size and simulation time,
the machine learning (ML) potential model is able to reproduce the experimental enthalpy profile increasing the hydration film and provides converged densities with a well-defined minimum between the first water layer and water molecules beyond. Further, time correlation function combined to ML method allows for the observation of dynamic occurrences on surfaces, including residence lifetime and the exchange of water molecules between hydration layers. In addition to offering valuable insights into the hydration mechanisms of kaolinite surfaces, for the first time we validate the use of ML force field to overcome the limitations of both first principles and classical molecular dynamics simulations for investigating adsorption phenomena.

January 13, 2025

# I. Introduction

Kaolinite ($Al_2Si_2O_5(OH)_4$) is a common clay mineral formed from altered rocks and a source of alumina impurities.¹,² Known for its excellent adsorbent properties, it has been largely studied for the removal of pollutants, such as heavy metals,³,⁴ dyes⁵ and phenols⁶,⁷ from wastewaters, thus constituting a promising material for the environmental domain. Yet, due to its main composing elements (Al and Si), kaolinite is also known for being deleterious to metallurgical processes,⁸ reducing both efficiency and productivity,⁹,¹⁰ especially in iron ore recovery. Hence, the impurities contained in minerals such as kaolinite together with other silicates - i.e., phosphates or sulfides - have to be removed from ores by rejecting kaolinite.¹¹ Of note, the removal of kaolinite from iron ores is challenging.¹¹,¹² Although magnetic and gravity separation are the common beneficiation routes of iron ores, the decrease in liberation size limits the performance of the beneficiation process. Nowadays, the most effective method to upgrade iron ore concentrates is the reverse flotation,¹³ a separation technique based on the selective adsorption of reagents onto the mineral surfaces to tune the surface tension of the minerals.

Both pollutant removal using kaolinite and beneficiation processes to remove it share the
same hydrated conditions during the performances of the process. Indeed, in all mineral
processing routes, particularly in flotation, ores are mixed with water for transport, for dust
management, and for most of the processing stages. Thus, prior to any pollutant capture by
adsorption or surface tension change by adsorption of flotation reagents, understanding the
interactions between the mineral surfaces and water is of paramount interest.¹⁴

Recently, several studies have focused on the hydration mechanism of kaolinite surfaces. ¹⁵⁻¹⁸
Molecular simulations conducted on the basal hydroxyl surface have established that a sole
water molecule adsorbs by forming three hydrogen bonds with the surface. ¹⁹ The first two
bonds involve the oxygen atoms of water and the surface protons, while the third one is
established between the water hydrogen atom and the oxygen atom of the hydroxyl termi-
nation on surface. When one increases the coverage, water molecules generate a monolayer
held together by hydrogen bonds and described as an ice-like structure. ²⁰ Furthermore, the
wettability of the kaolinite surfaces depends upon the amphoteric properties of the hydrox-
ylated surfaces, as well as its physicochemical properties and roughness. ²¹ Albeit methods
such as the sessile drop technique are able to assess the water contact angle,²² the exper-
imental measure of adsorption energy of water at interfaces is very challenging whenever
atomic resolution is desired. Computer simulations, especially when extended, provide ac-
cess to detailed mechanistic information on solvent configurations and may ascertain crucial
dynamical events over the adsorption process.

The goal of the present work is a comprehensive investigation of the adsorption mechanisms
of water onto both the tetrahedral (siloxane) and the octahedral (hydroxyl) basal surfaces of
kaolinite.²³,²⁴ Theoretical modeling relies first of all on classical molecular dynamics (CMD)
and ab initio molecular dynamics (AIMD) simulations. The former requires a reasonable
time scale to satisfy the ergodic hypothesis,²⁵ yet it may overcome the limits of infinitesimal
time and space during the adsorption since based on the principles of Newtonian mechanics.
The latter explicitly evaluates the full electronic structure using density functional the-
ory (DFT),²⁶ effectively gauging the reactive events at the interface. However, despite the
progress made for both classical<sup>27</sup> and ab initio simulations,<sup>28</sup> when it comes to interpreting
the structural properties of water/solid interactions provided by the experiments, CMD fail
quite often in providing the microscopic interface picture.<sup>29,30</sup> On the other hand, although
first-principles techniques would allow for extended simulations, their high computational
cost represents a severe bottleneck in gauging complex interplay of many phenomena at the
interface, such as surface diffusion or water exchange. Thus, AIMD simulations are limited to
short times of the order of a few picoseconds and are inherently restricted to relatively small
system sizes.<sup>31</sup> First-principles-based potential energy surfaces<sup>32</sup> may extend the simulation
time to nanoseconds and have been demonstrated to be quite effective in predicting adsorp-
tion energies,<sup>33</sup> especially when active learning techniques to generate the training data set
are used.<sup>34-37</sup> Although nowadays approaches capable of automatically selecting configura-
tions in many different applications exist,<sup>34</sup> strategies to lighten the training cost are still in
their infancy. Here, a machine learning force field (MLFF)<sup>38,39</sup> model was iteratively trained
based on AIMD trajectories to generating new data points in CMD trajectories. The hy-
drated MLFF model was employed to evaluate the isosteric enthalpies of water molecules
on kaolinite surfaces at large surface coverages and contrasted to the experimental profile.
Computational details are presented in the next section, while the following sections provide
the results, the associated discussion, and the summarizes of the main findings, including
the validation of the MLFF model.

II. Theoretical Modeling

Classical simulations

CMD was performed using the open Large-scale Atomic Molecular Massively Parallel Sim-
ulator (LAMMPS) package<sup>40</sup> in the NVT ensemble (constant particle number, volume, and
300 K temperature over 100 ns). The equations of motion were solved using the velocity-
Verlet integrator with a timestep of 1 fs, under the Nosé-Hoover thermostat.<sup>41,42</sup> A 3×2×1
triclinic supercell containing 204 atoms was generated. To prevent unwanted interaction be-

<table><caption>Table 1: Parameters of water atoms in the SPC/E model.</caption><thead><tr><th colspan="4">Non-bonding Lennard-Jones potential</th></tr><tr><th>Species</th><th>q(e)</th><th>ε(kcal · mol<sup>-1</sup>)</th><th>σ(Å)</th></tr></thead><tbody><tr><td>O</td><td>-0.8476</td><td>0.1553</td><td>3.166</td></tr><tr><td>H</td><td>0.4238</td><td>0.0000</td><td>0.0000</td></tr><tr><th colspan="4">Bond-Harmonic potential</th></tr><tr><td>Bond</td><td>r<sub>0</sub> (Å)</td><td>K(kcal · mol<sup>-1</sup>)</td><td></td></tr><tr><td>OH</td><td>1.0</td><td>554.13</td><td></td></tr><tr><th colspan="4">Angle - Harmonic potential</th></tr><tr><td>Angle</td><td>K(kcal · mol<sup>-1</sup> · Rad<sup>-2</sup>)</td><td>θ<sub>0</sub>(Å)</td><td></td></tr><tr><td>H-O-H</td><td>45.77</td><td>109.47</td><td></td></tr></tbody></table>

tween the water molecules and the reciprocal surface, a vacuum of 30 Å was added on top of the interface. Then, the relative volume was filled with water molecules using the Packmol software.43,44 Table 1 lists the parameters dealing with the 3-site rigid SPC/E model45 used to characterize water molecules. Table 2 lists the parameters of the modified Clay Force Field which contains an explicit term to accurately describe the Al-O-H angle as well as the OH bonds.16,46

<table><caption>Table 2: Parameters of the modified Clay Force Field used for CMD simulations.</caption><thead><tr><th colspan="4">Non-bonding Lennard-Jones potential</th></tr><tr><th>Species</th><th>q(e)</th><th>ε(kcal · mol<sup>-1</sup>)</th><th>σ(Å)</th></tr></thead><tbody><tr><td>Tetrahedral Si</td><td>2.100</td><td>1.8405 × 10<sup>-6</sup></td><td>3.3020</td></tr><tr><td>Octahedral Al</td><td>1.5750</td><td>1.3298 × 10<sup>-6</sup></td><td>4.2718</td></tr><tr><td>Bridging O</td><td>-1.0500</td><td>0.1554</td><td>3.1655</td></tr><tr><td>Hydroxyl O</td><td>-0.9500</td><td>0.1554</td><td>3.1655</td></tr><tr><td>Hydroxyl H</td><td>0.4250</td><td>0.000</td><td>0.0000</td></tr><tr><th colspan="4">Bonds - Morse potential</th></tr><tr><td>Bond</td><td>D<sub>0</sub>(kcal · mol<sup>-1</sup>)</td><td>A(Å<sup>-1</sup>)</td><td>R<sub>0</sub>(Å)</td></tr><tr><td>OH</td><td>2.1000</td><td>1.8405 × 10<sup>-6</sup></td><td>3.3020</td></tr><tr><th colspan="4">Angle - Harmonic potential</th></tr><tr><td>Angle</td><td>K(kcal · mol<sup>-1</sup> · Rad<sup>-2</sup>)</td><td>θ<sub>0</sub>(Å)</td><td></td></tr><tr><td>Al-O-H</td><td>15</td><td>110</td><td></td></tr></tbody></table>
## Ab initio simulations

DFT calculations at 0 K and AIMD simulations at 300 K were performed on a 3×1×1 kaolin-ite cell. A vacuum of 10 Å was systematically maintained between the interface and the upper limit of the cell. Further, to prevent any drift of the kaolinite slab during the simulations, the atoms corresponding to the first layer of the cell were frozen to their equilibrium position. Total energy and structure of the system were determined using the Vienna Ab initio Simulation Program (VASP).<sup>47</sup> The semi-local Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional in the generalized gradient approximation (GGA)<sup>48</sup> was employed and the projector augmented wave (PAW)<sup>49,50</sup> was used to describe the electron-ion interactions. The Kohn-Sham equations<sup>51</sup> were self-consistently<sup>52</sup> solved until a threshold energy difference of 10<sup>-7</sup> eV between two consecutive cycles was reached. The plane wave cut-off energy showed convergence and set at 500 eV. To further the total energy convergence, a Methfessel Paxton<sup>53</sup> smearing of σ=0.1 eV was applied. The structural relaxation was performed until the forces were within 0.01 eV/Å. The Brillouin Zone integration was performed using only Γ-point because of the cell dimension. Van der Waals forces were integrated over the calculations using the D3 correction of Grimme.<sup>54</sup> The AIMD simulations in the NVT ensemble (300 K over 100 ps, timestep of 1 fs) were performed with an energy cut-off of 400 eV, a threshold energy difference between two cycles of 10<sup>-5</sup> eV and a thermalization period of 10 ps. Structural analysis for a single water molecule adsorbed upon the slab was conducted by performing an AIMD simulation for 200 ps. The system's final arrangement was subsequently optimized at 0 K, with both the initial and optimized configurations for the slab and slab-water interface being available in Supporting Information (SI) for the purpose of reproducibility of the simulations.

## Machine learning simulations

MLFF as implemented in VASP 6.4<sup>38,39</sup> was initially trained using 2 × 10<sup>4</sup> AIMD steps. As shown in Figure S1, the values for both radial and angular cutoffs dealing with the
descriptor were estimated by reducing the Bayesian error for the forces. At this stage, the active learning procedure outlined by Schran et al.³⁶ was employed to progressively refine the MLFF model. After a preliminary training, the model was immediately used to run classical simulations in order to extrapolate and generate data points beyond the portion of potential energy surface already explored. Configurations with the largest energy variance during the prediction were systematically chosen to fill the unexplored configurational space and boost the data set. The training ended whenever the root mean square error (RMSE) between the calculated forces from MLFF and the explicit DFT was within a threshold. Following the completion of the training, the data collected during the production run were used to determine the observables whereas statistical Mann-Kendall tests⁵⁵ were performed over the equilibration period. As presented in figures S2-S5, the MLFF model was validated by comparing the profiles of the radial distribution function computed for the pairs of atoms to the same profiles obtained from AIMD calculations. This comparison was systematically performed starting from the slab up to the fully hydrated system. Table 3 outlines the comparison between AIMD and MLFF in terms of adsorption energy values obtained for one water molecule, as well as half and full hydrated system. Figure S6 shows the IR spectrum obtained through MLFF dynamics at 100 K. If contrasted to the experimental one,⁵⁶ the theoretical model successfully replicates the presence of the four OH stretching modes, albeit displaying a minor deviation in their positions. Finally, figure S7 illustrates the computational cost associated with the AIMD method in contrast to the trained MLFF potential across a variety of structures differing in dimensions. The total computational time needed by AIMD shows a significant increase with the size of the system. In contrast, the total time necessary to predict the electronic structure using MLFF is considerably decreased, highlighting a linear correlation with the size of the system. For additional analysis and discussion regarding the ML computational cost, the reader is directed to the SI document.
<table><caption>Table 3: Comparison between AIMD and MLFF dealing with the adsorption energy values in kJ·mol⁻¹.</caption><thead><tr><th></th><th>AIMD</th><th>MLFF</th></tr></thead><tbody><tr><td>Kaolinite + 1 H<sub>2</sub>O</td><td>-60.50</td><td>-58.54</td></tr><tr><td>Kaolinite + 1/2 monolayer H<sub>2</sub>O</td><td>-63.85</td><td>-64.05</td></tr><tr><td>Kaolinite + 1 monolayer H<sub>2</sub>O</td><td>-69.33</td><td>-68.03</td></tr></tbody></table>

## Evaluation of the properties

The adsorption energies were calculated according to the expression:

$$ \Delta E_{ads} = E_{SM} - E_{S} - E_{M} \qquad (1) $$

where $E_{SM}$, $E_{S}$, and $E_{M}$ refer to the potential energies of the relaxed molecule interacting with the surface, of the clean surface and of the molecule in the gas phase, respectively. Similarly, the calculations of the dispersion forces followed the expression below, in which the considered energies are dispersive energies:

$$ \Delta E^{d} = E_{SM}^{d} - E_{S}^{d} - E_{M}^{d} \qquad (2) $$

Besides, when several water molecules were considered, the adsorption energies per molecule were evaluated including the dispersion energies as follows:

$$ \Delta E = \frac{E_{SM} - E_{S} - n \times E_{M}}{n} \qquad (3) $$

where $n$ is the amount of adsorbed water molecules. The isosteric enthalpy dealing with the adsorption of $n$ molecules was estimated from the adsorption energies of the system before ($\Delta E_{ads}^{1}$) and after ($\Delta E_{ads}^{n+1}$) the adsorption of those $n$ molecules as follows:

$$ \Delta H = \frac{\Delta E_{ads}^{1+n} - \Delta E_{ads}^{1}}{n} \qquad (4) $$
Finally, for gauging the kinetics of water exchange between the first and second (or higher) water layers, MLFF was used with the stable-state-picture (SSP) technique.⁵⁷ Accordingly, the time correlation function was obtained as

$$C(t) = 1 - \frac{1}{N} \sum_{i}^{N} \langle s_{r}^{(i)}(0) s_{p}^{(i)}(t) \rangle \quad (5)$$

where the sum ran over all water molecules N settled in the first solvation layer at t = 0 and where $s_{r}^{(i)}$ and $s_{p}^{(i)}$ are step functions corresponding to reactant (product) state if the water molecule i was in the first (second or higher) solvation layer.

# III. Results and discussion

## Single water molecule

Figure 1 presents the positions tested for the adsorption of a single water molecule on the tetrahedral surface at 0 K. For each configuration, the water molecule was oriented with the hydrogen pointing towards or away from the surface. Among them, the most favored

<figure><img src="image_3.png" /><figcaption>Figure 1: Top view dealing with the positions tested for a single water molecule on the tetrahedral surface of kaolinite.</figcaption></figure>
configuration turned out to be the one with the water molecule at the center of a Si-O ring (Figure 2), with the water hydrogen atoms pointing towards the surface ($d = 2.2$ Å). The adsorption energy in this position was $-37.6 \text{ kJ} \cdot \text{mol}^{-1}$ with a dispersion contribution of $-13.8 \text{ kJ} \cdot \text{mol}^{-1}$. Next, three configurations were investigated in order to analyze the dissociation of

<figure><img src="image_4.png" /><figcaption>Figure 2: The most favoured configuration dealing with the adsorption of a single water molecule onto the tetrahedral surface of kaolinite, with bonds lengths (left) and the difference of charge density (right). Isosurfaces span from $10^{-3}$ electrons/Å$^3$ for the yellow area to $-10^{-3}$ electrons/Å$^3$ for the blue area.</figcaption></figure>

the water molecule on the surface. In the *close* configuration ($C_c$), both ions were separated at distances ranging from 1.6 to 3.2 Å on top of the same Si-O ring. In the *medium* configuration ($C_m$), the ions were set in two different Si-O rings at ~5 Å from each other. Finally, in the *far* configuration ($C_f$), both ions were separated by at least one Si-O ring at 9.5 Å from each other. The vertical distance from the surface was maintained over the simulations at 1.6 Å for both ions. Both $C_c$ and $C_m$ conducted to the reformation of the water molecule, with adsorption energies ranging from $-37.6$ to $-32.0 \text{ kJ} \cdot \text{mol}^{-1}$ and dispersion contributions ranging from $-13.5$ to $-10.7 \text{ kJ} \cdot \text{mol}^{-1}$. Albeit for $C_f$ the proton was adsorbed onto the surface whereas the hydroxyl was expelled, the associated enthalpy was found highly endothermic. Hence, water molecule does not dissociate when adsorbed onto the tetrahedral surface.

Then, we moved to study the octahedral interface. Figure 3 presents the relative basal (001) kaolinite surface used to perform the theoretical simulations. Consistently with previous cal-
<figure><img src="image_5.png" /><figcaption>Figure 3: Kaolinite (001) basal surface used for theoretical modeling and shown from different angles. Parameters of the triclinic cell were taken from Bish and co-workers⁵⁸ and are as follows: a = 5.15 Å, b = 8.94 Å and c = 7.39 Å; α = 91.93°, β = 105.05° and γ = 89.80°.</figcaption></figure>

culations,¹⁹ the water molecule adsorbed with an energy of -62.4 kJ·mol⁻¹ and a dispersion contribution of -17.7 kJ·mol⁻¹, values significantly more negative than for the tetrahedral surface. As shown in Figure 4, the water molecule formed three hydrogen bonds with the kaolinite octahedral surface: two bonds between the oxygen atom of the water molecule and the hydrogen atoms of the surface ( 1.95 < *d* < 2.1 Å) and one bond between a water hydrogen atom and a surface oxygen atom (1.7 < *d* < 1.73 Å). Of note, the latter distance values are consistent with the ones calculated by Tunega, Chen and their co-workers.¹⁵,¹⁹ When it came to investigating the dissociation of the water molecule on this surface, we found that the proton and the hydroxyl recombine only when their distance is within 3 Å.
<figure><img src="image_6.png" /><figcaption>Figure 4: Left: side and top view dealing with the adsorption of a single water molecule onto the octahedral surface of kaolinite. Right: difference of charge densities for of the water molecule adsorbed onto the surface. Isosurfaces span from $10^{-3}$ electrons/Å$^3$ for the yellow area to $-10^{-3}$ electrons/Å$^3$ for the blue area.</figcaption></figure>

## From single water molecule to the monolayer

Next, water molecules were continuously added onto both the tetrahedral and the octahedral surfaces, until a first layer of adsorbed molecules was reached (corresponding to an adsorbed amount of ~ 0.20 mmol·g$^{-1}$). As shown in Figure 5, on the octahedral surface, water formed a hexagonal-shaped network of bonded molecules located at approximately 3 Å from the surface. By contrast, on the tetrahedral surface, water molecules were maintained at the same axial distance of 3 Å and did not exhibit this hexagonal framework. Figure S8 presents the energy profiles obtained from DFT calculations at 0 K onto both surfaces. Overall, water was more favourably adsorbed onto the octahedral surface, which was consistent with our findings for a single water molecule but also consistent with the rather hydrophilic behaviour of this surface. The enthalpy of adsorption decreased from -62.4 kJ·mol$^{-1}$ for a single water molecule to -86.9 kJ·mol$^{-1}$ for the monolayer, which meant that the water molecules were stabilised by the formation of a network on the surface. Beyond the formation
of the monolayer, corresponding to the minimum in the enthalpy profile, the energy increased again with coverage and at higher values (> 0.6 mmol/g), and stabilized around -70 kJ·mol⁻¹. Consistently with what was previously observed by Tunega et al.,¹⁹ a smaller enthalpy value of -65.5 kJ·mol⁻¹ was found for the monolayer formed onto the tetrahedral surface. When more water molecules were added on top of this surface, they condensed in the gas phase without forming any structural arrangement. On the octahedral surface, a second water

<figure><img src="image_7.png" /><figcaption>Figure 5: Monolayer on top of kaolinite's octahedral surface seen from (a) the b-axis (side view), from (b) the c-axis (top view) and from (c) the c-axis with the monolayer represented by dotted lines.</figcaption></figure>
layering formed at approximatively 6 Å from the surface. Beyond that, water molecules
formed bulk-like water, which raised the question of the existence of the third layer.¹⁵

<figure><img src="image_8.png" /></figure>

Figure 6: (a) Isosteric enthalpy of the adsorption of water onto kaolinite octahedral surface estimated from CMD, AIMD and MLFF simulations at 300 K. (b) Increasing the adsorbed amount of water beyond the monolayer, comparison between the adsorption energy profiles obtained from the experimental dual-site Langmuir-Freundlich (DSLF) model (blue solid curve<sup>59</sup>) and the MLFF model (black filled circles). (c) Comparison between AIMD and MLFF techniques dealing with oxygen water densities. By convention, the interfacial layer is settled to 0 Å.

Beyond the limits of CMD and AIMD: MLFF

To address this topic, we had to reasonable scale both size and simulation time in order to
improve the statistical sampling. Since the adsorption of water does not involve a chemical
reaction, we first compared at room temperature the quality of CMD to AIMD simulations.
Figure 6(a) illustrates the enthalpy profiles derived on the octahedral surface through CMD,
AIMD, and MLFF simulations. Overall, the energy differences of ~ 10 kJ·mol⁻¹ could be
attributed to the looser network formed by water molecules during the CMD trajectory. The
MLFF model successfully replicates the AIMD profile, albeit with a slight shift in the curve’s
minimum and an underestimation of the maximum value. Additionally, as the MLFF profile
covers a broader range, it converges to a lower energy level by 2 kJ·mol⁻¹ compared to the
AIMD profile. Besides, Figure S9 shows a pictorial representation increasing the amount
of adsorbed water molecules over the CMD simulation. The minimum energy of the CMD
profile corresponded to -62.0 kJ·mol⁻¹ and an adsorbed amount of ~ 0.15 mmol·g⁻¹. At
higher values of surface coverage, the enthalpy stabilized at -47 kJ·mol⁻¹. With regard
to AIMD simulations, Figure 7 displays a few snapshots obtained by progressively adding
water molecules onto the octahedral surface. The monolayer in this case corresponded to an
adsorbed amount comprised between 0.15 and 0.20 mmol·g⁻¹, with an adsorption energy of
-73.6 kJ·mol⁻¹ and located at 3 Å from the surface. Beyond the monolayer, the enthalpy

<figure><img src="image_9.png" /><figcaption>Figure 7: Snapshots dealing with the adsorption of water onto kaolinite octahedral surface at 300 K over AIMD simulations. Values are in mmol/g.</figcaption></figure>
calculated by AIMD simulations increased with coverage and stabilized at -53.4 kJ·mol⁻¹.
As shown in Figure 6b, this value was consistent with the adsorption isotherms experiments
for low adsorbed amounts⁵⁹ and comparable with the adsorption/desorption energies com-
monly observed on metals and minerals surfaces.⁶⁰,⁶¹ To improve the simulation accuracy for

<figure><img src="image_10.png" /><figcaption>Figure 8: Time correlation functions obtained from MLFF simulations describing the resi-
dence lifetime of water molecules in the first solvation layer at the kaolinite surface. For the
sake of pictorial representation, dimensions of Si (blue) and Al (gray) have been magnified in
comparison to O (red) and H (white) to accentuate the distinct stepped-edge configuration
of the computational models.</figcaption></figure>

larger amounts of water, a first-principles-trained potential energy surface was obtained and
compared to the experimental isotherm.⁵⁹ While showing a smoother profile (Figure 6b),
the MLFF model was able to reproduce the experimental one, reaching the plateau for high
values of surface coverages at ~ -39 kJ·mol⁻¹, nearby the enthalpy of condensation of water
(-40.7 kJ/mol). Figure 6c shows how further insights on water structuring were obtained by
monitoring the water oxygen atoms density above the kaolinite slab, settled at 0 Å. First, the
ML simulations qualitatively recovered the explicit profile obtained by AIMD simulations,
exhibiting a pronounced peak corresponding to the first adsorption layer characterized by a
highly structured water arrangement. Noteworthy, the AIMD simulations failed in providing well defined layering structures at higher distances from the slab, most probably due to the poor statistical sample. By increasing the order of magnitude by a factor three, the data collected by the trajectory were able to design a well converged density profile, exhibiting water structuring in the first two layers. The first minimum in the density profile indicated the poor water exchange within the first two hydration films whereas the second minimum precedes the formation of bulk water, rather complete beyond 8 Å.

Finally, an investigation was conducted to assess the predictive capability of the MLFF model in anticipating the damping effect on the surfactant water motion caused by the presence of edge surfaces, which are accountable for the amphoteric nature of kaolinite.⁶² Within the confines of this theoretical framework, the AIMD simulations encountered limitations in terms of size and duration whereas CMD simulations failed to offer essential insights at the edge interface due to its notably high reactivity. Particularly in a hydrated environment, silanol and aluminol groups on the edges exhibit hydroxyl groups coordinated to only one silicon or aluminium atom.⁶³ Consequently, the edge surface is distinguished by a considerable quantity of exceedingly reactive terminations that have the potential to attenuate water motion upon contact with the surface, thereby influencing the overall interfacial reactivity of kaolinite.⁶² Starting from the octahedral basal surface delineated in Figure 8(a), two models were devised to incorporate and gradually augment the edge surface area. The models showcased in Figures 8(b) and 8(c) corresponded to edge coverages of 36% and 43% relative to the basal slab, respectively. By employing the time correlation function $C(t)$ in conjunction with MLFF, insights into the kinetics and exchange of water molecules between the first and second layers were provided. The correlation profiles illustrated in Figure 8 indicate that the basal model demonstrates the briefest residence lifetime of water molecules within the primary hydration film. Conversely, it is evident from models (b) and (c) that the residence lifetime was influenced by the edge interface and its degree of hydroxylation. Notably, for high hydroxylation values, -Si-OH and -AlOH₂OH terminations exhibit amphoteric charac-
teristics and high reactivity.⁶⁴ Subsequently, the terminations may transition gradually from
–OH₂ to –O via –OH by releasing their protons, contingent on the pH of the solution.⁶⁵,⁶⁶
On the whole, correlation profiles at heightened hydroxylation values (≡ R–OH) exhibited a
shorter residence lifetime within the first layer compared to lower values (≡ R–O), suggesting
that the dangling terminations on the edge may also impact the water molecules’ kinetics at
the slab interface.

IV. Concluding Remarks

In this study, a comprehensive theoretical analysis of the hydration mechanisms occurring on
the surfaces of kaolinite has been presented. First-principles simulations allowed identifying
the predominant structural water organizations involved at the interface, revealing a highly
organized water arrangement within the initial hydration layer. The employment of machine
learning extended simulations facilitated the exploration of water density profiles and the
observation of dynamic processes, such as water exchange between the first hydration layer
and subsequent water layers. Findings in this direction could offer valuable insights into the
structure-property relationship at these interfaces and pave the way for more accurate ad-
sorption models and improved characterization methods within the flotation industry. Yet,
due to the intricate interplay of various phenomena at the solid-liquid interface, including
surface diffusion and fluctuations in hydrogen bonding, understanding the origins and mech-
anisms of specific microscopic properties remains a significant challenge. Studies aimed at
assessing the role of selective adsorption of collectors and depressants as well as at improving
the transfer of ML model strategies in the contest of mineral/water interfaces are currently
underway in our group.
Acknowledgement

Part of this work has been conducted in the frame of the industrial chair MULTIMINE
funded by ArcelorMittal, Région Grand Est, Metz Métropôle and Université de Lorraine.
This work was granted access to the HPC resources of TGCC under the allocation 2023-
A0140810433 made by GENCI. We also thank Professor Tomáš Bučko from Comenius Uni-
versity, Bratislava for support and discussions.

Supporting Information Available

Further details on both calculation of the isosteric enthalpy of adsorption from water isotherms and validation of the MLFF potential are available in the supplementary material. Additional data are available from the corresponding authors upon reasonable request.

References

(1) Johnston, C. J.; Pepper, R. A.; Martens, W. N.; Couperthwaite, S. Improvement of aluminium extraction from low-grade kaolinite by iron oxide impurities: Role of clay chemistry and morphology. *Miner. Eng.* **2022**, *176*, 107346.

(2) Nnanwube, I. A.; Keke, M.; Onukwuli, O. D. Assessment of Owhe kaolinite as poten-
tial aluminium source in hydrochloric acid and hydrogen peroxide solutions: Kinetics
modeling and optimization. *Clean. Chem. Eng.* **2022**, *2*, 100022.

(3) Gupta, S. S.; Bhattacharyya, K. G. Adsorption of heavy metals on kaolinite and mont-
morillonite: a review. *Phys. Chem. Chem. Phys.* **2012**, *14*, 6698.

(4) Mustapha, S.; Ndamitso, M.; Abdulkareem, A.; Tijani, J.; Mohammed, A.; Shuaib, D. Potential of using kaolin as a natural adsorbent for the removal of pollutants from tannery wastewater. *Heliyon* **2019**, *5*, e02923.
(5) Khan, T. A.; Khan, E. A.; Shahjahan, Removal of basic dyes from aqueous solution by adsorption onto binary iron-manganese oxide coated kaolinite: Non-linear isotherm and kinetics modeling. *Appl. Clay Sci.* **2015**, *107*, 70–77.

(6) Alkaram, U. F.; Mukhlis, A. A.; Al-Dujaili, A. H. The removal of phenol from aqueous solutions by adsorption using surfactant-modified bentonite and kaolinite. *J. Hazard. Mater.* **2009**, *169*, 324–332.

(7) Abdel-Aleem, A. A. M.; Abdel-Tawab, M. S.; Hassouna, M. E.-K. M. Phenol removal from aqueous environments by natural &amp chemically modified kaolin clay. *Environ. Qual. Manag.* **2022**, *32*, 119–135.

(8) Rivera, O.; Pavez, O.; Kao, J. L.; Nazer, A. Metallurgical characterization of kaolin from Atacama, Chile. *REM - Int. Eng. J.* **2016**, *69*, 473–478.

(9) Zhan, W.; Liu, Y.; Shao, T.; Han, X.; Pang, Q.; Zhang, J.; He, Z. Evaluating the Effect of MgO/Al2O3 Ratio on Thermal Behaviors and Structures of Blast Furnace Slag with Low Carbon Consumption. *Crystals* **2021**, *11*, 1386.

(10) Agrawal, A.; Das, K.; Singh, B. K.; Singh, R. S.; Tripathi, V. R.; Kundu, S.; Padmapal, R.; Singh, M. K. Means to cope with the higher alumina burden in the blast furnace. *Iron & Steel.* **2020**, *47*, 238–245.

(11) Rodrigues, O. M. S.; Peres, A. E. C.; Martins, A. H.; Pereira, C. A. Kaolinite and hematite flotation separation using etheramine and ammonium quaternary salts. *Miner. Eng.* **2013**, *40*, 12–15.

(12) Fan, G.; Wang, L.; Cao, Y.; Li, C. Collecting Agent–Mineral Interactions in the Reverse Flotation of Iron Ore: A Brief Review. *Minerals* **2020**, *10*, 681.

(13) Filippov, L.; Severov, V.; Filippova, I. An overview of the beneficiation of iron ores via reverse cationic flotation. *Int. J. Miner. Process.* **2014**, *127*, 62–69.
(14) Foucaud, Y.; Lebègue, S.; Filippov, L. O.; Filippova, I. V.; Badawi, M. Molecular Insight into Fatty Acid Adsorption on Bare and Hydrated (111) Fluorite Surface. *J. Phys. Chem. B* **2018**, *122*, 12403–12410, Number: 51.

(15) Chen, J.; Min, F.-f.; Liu, L.-y.; Liu, C.-f. Mechanism research on surface hydration of kaolinite, insights from DFT and MD simulations. *Appl. Surf. Sci.* **2019**, *476*, 6–15.

(16) Pouvreau, M.; Greathouse, J. A.; Cygan, R. T.; Kalinichev, A. G. Structure of Hydrated Kaolinite Edge Surfaces: DFT Results and Further Development of the ClayFF Classical Force Field with Metal–O–H Angle Bending Terms. *J. Phys. Chem. C* **2019**, *123*, 11628–11638, Number: 18.

(17) Yavna, V.; Nazdracheva, T.; Morozov, A.; Ermolov, Y.; Kochur, A. Ab Initio Simulation of the IR Spectrum of Hydrated Kaolinite. *Crystals* **2021**, *11*, 1146.

(18) Miao, Y.; Yan, H.; Qiu, X.; Zhou, X.; Zhu, D.; Li, X.; Qiu, T. Adsorption of hydrated Al3 on the kaolinite (001) surface: A density functional theory study. *Appl. Clay Sci.* **2022**, *223*, 106498.

(19) Tunega, D.; Gerzabek, M. H.; Lischka, H. Ab Initio Molecular Dynamics Study of a Monomolecular Water Layer on Octahedral and Tetrahedral Kaolinite Surfaces. *J. Phys. Chem. B* **2004**, *108*, 5930–5936, Number: 19.

(20) Hu, X. L.; Michaelides, A. Water on the hydroxylated (001) surface of kaolinite: From monomer adsorption to a flat 2D wetting layer. *Surf. Sci.* **2008**, *602*, 960–974, Number: 4.

(21) Hu, X. L.; Michaelides, A. Ice formation on kaolinite: Lattice match or amphotерism? *Surf. Sci.* **2007**, *601*, 5378–5381, Number: 23.

(22) Shang, J.; Flury, M.; Harsh, J. B.; Zollars, R. L. Contact angles of aluminosilicate
clays as affected by relative humidity and exchangeable cations. *Colloids Surf. A: Physicochem. Eng. Asp.* **2010**, *353*, 1–9, Number: 1.

(23) Bleam, W. F. Atomic theories of phyllosilicates: Quantum chemistry, statistical mechanics, electrostatic theory, and crystal chemistry. *Rev. Geophys.* **1993**, *31*, 51–73.

(24) Murray, H. H. Traditional and new applications for kaolin, smectite, and palygorskite: a general overview. *Appl. Clay Sci.* **2000**, *17*, 207–221.

(25) Coveney, P. V.; Wan, S. On the calculation of equilibrium thermodynamic properties from molecular dynamics. *Phys. Chem. Chem. Phys.* **2016**, *18*, 30236–30240.

(26) Marx, D.; Hutter, J. *Ab Initio Molecular Dynamics*; Cambridge University Press, 2009.

(27) Cisneros, G. A.; Wikfeldt, K. T.; Ojamäe, L.; Lu, J.; Xu, Y.; Torabifard, H.; Bartók, A. P.; Csányi, G.; Molinero, V.; Paesani, F. Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential Energy Functions. *Chem. Rev.* **2016**, *116*, 7501–7528.

(28) Gillan, M. J.; Alfè, D.; Michaelides, A. Perspective: How good is DFT for water? *J. Chem. Phys.* **2016**, *144*, 130901.

(29) Quaranta, V.; Hellström, M.; Behler, J. Proton-Transfer Mechanisms at the Water–ZnO Interface: The Role of Presolvation. *J. Phys. Chem. Lett.* **2017**, *8*, 1476–1483.

(30) Jakub, Z.; Meier, M.; Kraushofer, F.; Balajka, J.; Pavelec, J.; Schmid, M.; Franchini, C.; Diebold, U.; Parkinson, G. S. Rapid oxygen exchange between hematite and water vapor. *Nat. Commun.* **2021**, *12*, 6488.

(31) He, X.; Zhu, Y.; Epstein, A.; Mo, Y. Statistical variances of diffusional properties from ab initio molecular dynamics simulations. *npj Comput. Mater.* **2018**, *4*, 18.

(32) Behler, J. Perspective: Machine learning potentials for atomistic simulations. *J Chem. Phys.* **2016**, *145*, 170901.
(33) Ras, E.-J.; Louwerse, M. J.; Mittelmeijer-Hazeleger, M. C.; Rothenberg, G. Predicting adsorption on metals: simple yet effective descriptors for surface catalysis. *Phys. Chem. Chem. Phys.* **2013**, *15*, 4436–4443.

(34) Smith, J. S.; Nebgen, B.; Lubbers, N.; Isayev, O.; Roitberg, A. E. Less is more: Sampling chemical space with active learning. *J. Chem. Phys.* **2018**, *148*, 241733.

(35) Bernstein, N.; Csányi, G.; Deringer, V. L. De novo exploration and self-guided learning of potential-energy surfaces. *npj Computational Materials* **2019**, *5*, 99.

(36) Schran, C.; Thiemann, F. L.; Rowe, P.; Müller, E. A.; Marsalek, O.; Michaelides, A. Machine learning potentials for complex aqueous systems made simple. *Proc. Natl. Acad. Sci.* **2021**, *118*, e2110077118.

(37) Vandermause, J.; Xie, Y.; Lim, J. S.; Owen, C. J.; Kozinsky, B. Active learning of reactive Bayesian force fields applied to heterogeneous catalysis dynamics of H/Pt. *Nat. Commun.* **2022**, *13*, 5183.

(38) Ryosuke, J.; Ferenc, K.; Georg, K. On-the-fly machine learning force field generation: Application to melting points. *Phys. Rev. B* **2019**, *100*, 014105.

(39) Ryosuke, J.; Jonathan, L.; Ferenc, K.; Kresse, G.; Menno, B. Phase Transitions of Hybrid Perovskites Simulated by Machine-Learning Force Fields Trained on the Fly with Bayesian Inference. *Phys. Rev. Lett.* **2019**, *122*, 225701.

(40) Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; In 't Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comp. Phys. Commun.* **2022**, *271*, 108171.

(41) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. *Phys. Rev. A* **1985**, *31*, 1695–1697, Number: 3 Publisher: American Physical Society.
(42) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511–519, Number: 1 Publisher: American Institute of Physics.

(43) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157–2164, Number: 13.

(44) Martínez, J. M.; Martínez, L. Packing optimization for automated generation of complex system's initial configurations for molecular dynamics and docking. J. Comput. Chem. 2003, 24, 819–825, Number: 7.

(45) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269–6271, Number: 24.

(46) Greathouse, J. A.; Durkin, J. S.; Larentzos, J. P.; Cygan, R. T. Implementation of a Morse potential to model hydroxyl behavior in phyllosilicates. J Chem. Phys. 2009, 130, 134713, Number: 13.

(47) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 1993, 47, 558–561, Number: 1.

(48) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868, Number: 18.

(49) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953–17979, Number: 24.

(50) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999, 59, 1758–1775, Number: 3.

(51) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138, Number: 4A.
(52) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B **1996**, *54*, 11169–11186, Number: 16.

(53) Methfessel, M.; Paxton, A. T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B **1989**, *40*, 3616–3621, Number: 6 Publisher: American Physical Society.

(54) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. **2010**, *132*, 154104.

(55) Schiferl, S. K.; Wallace, D. C. Statistical errors in molecular dynamics averages. J. Chem. Phys. **1985**, *83*, 5203–5209.

(56) Balan, E.; Delattre, S.; Guillaumet, M.; Salje, E. K. Low-temperature infrared spectroscopic study of OH-stretching modes in kaolinite and dickite. Am. Min. **2010**, *95*, 1257–1266.

(57) Laage, D.; Hynes, J. T. On the Residence Time for Water in a Solute Hydration Shell: Application to Aqueous Halide Solutions. J. Phys. Chem B **2008**, *112*, 7697–7701.

(58) Bish, D. L. Rietveld Refinement of the Kaolinite Structure at 1.5 K. Clays and Clay Minerals **1993**, *41*, 738–744, Number: 6.

(59) Lainé, J. Multiscale approaches for the reverse flotation of iron ores: the case of kaolinite. PhD thesis, Université de Lorraine, Nancy, 2023.

(60) Thiel, P. A.; Madey, T. E. The interaction of water with solid surfaces: Fundamental aspects. Surf. Sci. Rep. **1987**, *7*, 211–385, Number: 6.
(61) Foucaud, Y.; Badawi, M.; Filippov, L. O.; Filippova, I. V.; Lebègue, S. Surface Properties of Fluorite in Presence of Water: An Atomistic Investigation. J. Phys. Chem. B **2018**, 122, 6829–6836, Number: 26.

(62) Liu, X.; Lu, X.; Wang, R.; Meijer, E. J.; Zhou, H.; He, H. Atomic scale structures of interfaces between kaolinite edges and water. Geochim. Cosmochim. Acta **2012**, 92, 233–242.

(63) Brady, P. V.; Cygan, R. T.; Nagy, K. L. Molecular Controls on Kaolinite Surface Charge. J. Colloid Interface Sci. **1996**, 183, 356–364, Number: 2.

(64) Schoonheydt, R.; Johnston, C. Dev. Clay Sci.; Elsevier, 2013; Vol. 5; pp 139–172.

(65) Liu, X.; Lu, X.; Sprik, M.; Cheng, J.; Meijer, E. J.; Wang, R. Acidity of edge surface sites of montmorillonite and kaolinite. Geochim. Cosmochim. Acta **2013**, 117, 180–190.

(66) Underwood, T.; Erastova, V.; Greenwell, H. C. Wetting Effects and Molecular Adsorption at Hydrated Kaolinite Clay Mineral Surfaces. J. Phys. Chem. C **2016**, 120, 11433–11449, Number: 21.
TOC Graphic

<figure><img src="image_11.png" /></figure>