# Investigation of Data Set Portability on Various Machine-Learned Interaction Potentials for Pyrophyllite Clay

Chloe Sanz, Colin Bousige, Pierre Mignon

► To cite this version:

Chloe Sanz, Colin Bousige, Pierre Mignon. Investigation of Data Set Portability on Various Machine-Learned Interaction Potentials for Pyrophyllite Clay. Journal of Physical Chemistry A, 2025, 130 (2), pp.581-588. 10.1021/acs.jpca.5c07319. hal-05534893

HAL Id: hal-05534893

https://hal.science/hal-05534893v1

Submitted on 3 Mar 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Investigation of Dataset Portability on Various Machine Learned Interaction Potentials for Pyrophyllite Clay

Chloe Sanz,† Colin Bousige,‡ and Pierre Mignon*,†

†Université Claude Bernard Lyon 1, CNRS, iLM UMR 5306, Villeurbanne, F-69100, France

‡Université Claude Bernard Lyon 1, CNRS, LMI UMR 5615, Villeurbanne, F-69100, France

E-mail: pierre.mignon@univ-lyon1.fr

## Abstract

The influence of the architecture and descriptors on the accuracy of a neural network potential is investigated. Particularly, we focused on high dimension machine learned interaction potentials (MLIPs) using descriptors based on atom centered symmetry functions, comprising embedding, attention masks or message passing procedures allowing to include information about the atomic neighbor environment. Four potentials trained on an identical representative dataset are validated, and their performances are compared against reference PBE-D3 results. All our MLIPs, regardless of the descriptors, are able to closely reproduce PBE-D3 energy and forces. They all have similar accuracy for the different parameters investigated (structural parameters, exfoliation energy, and vibrational spectra). This shows the importance of constructing a representative dataset in reaching the desired accuracy. In another effort to assess the accuracy, Grimme's D3 dispersion correction is added a posteriori to a MLIP fitted
on pure PBE data. This way, the dispersion correction is not affected by the cut-off
of the MLIP. This method produces accurate results for elastic parameters and vibra-
tional spectra, showing better accuracy for exfoliation energy and great stability for
structural parameters over various temperatures and cell sizes.

# 1 Introduction

The use of clay minerals permeates modern society. Their large surface area and layered
structure give them unique physical and chemical properties such as ionic exchange, swelling,
and adsorption properties.¹ Thus, they are used in various applications such as catalysis,²⁻⁴
play a role in the fabrication of plastics, paints, papers, or cosmetics,⁵ and can also be used
as depolluting agents.⁶,⁷ Pyrophyllite (Al₂(Si₄O₁₀)(OH)₂), is one of the simplest dioctahe-
dral clay (two-thirds of the octahedral sites occupied) with two SiO₄ tetrahedral (T) sheets
sandwiching one Al₂O₃ octahedral (O) sheet (see Fig. 1). The different TOT sheets are co-
valently bonded forming one layer, while the layers are held together through van der Waals
bonding. Contrary to other clay minerals, pyrophyllite does not have any substitution in the
T or O sheets and thus does not carry a permanent charge. Due to their great adsorption
capacity, pyrophyllite clays can have applications in the treatment of pollutants, e.g. for
removal of organic compounds like methylene blue⁸ and more recently for removal of nitrate
from soil and artificial fertilizers.⁹

To study pyrophyllite clay properties, in addition to experimental studies, theoretical
approaches such as atomistic simulations allow to explore a system's potential energy sur-
face (PES) allowing insight into structural, electronic, and dynamical properties.¹⁰ However,
atomistic simulations are confronted with a dilemma between highly accurate but compu-
tationally costly electronic structure methods, e.g. density functional theory (DFT), or
significantly more efficient but less accurate methods like classical potentials whose accu-
racy is force-field dependent. New approaches are now emerging, using machine learning
to approximate potential energy surfaces, allowing to run large-scale molecular dynamics
(MD) simulations with DFT-level accuracy for a fraction of the cost. These potentials are referred to as machine learned interaction potentials (MLIPs), and they are fitted from a dataset of reference data built from first-principle calculations.¹¹ The efficiency of the MLIP and the computational cost are influenced by i) its architecture and ii) by the descriptors (geometric inputs) used to encode the local environment, as such they define the dimensionality of the configuration space. A large choice of descriptors can be found in the literature, oftentimes alongside specific machine learning algorithms: atom-centered symmetry functions (ACSF),¹² Smooth Overlap of Atomic Positions,¹³ Atomic Cluster Expansion,¹⁴ etc. Many are predefined before training: they are often based on the local atomic environment, and usually include a cut-off truncating the interactions. For example, one can find ACSFs used in n2p2 package¹²,¹⁵ made of basis sets of Gaussian functions describing radial and angular environments, or the DeepPot-SE¹⁶ smooth descriptors including an additional embedding layer. Recently, the surge of graph neural networks and message passing methods has allowed to use learnable descriptors bypassing the need for explicit definition.¹⁷ Message passing allows to spread information from neighbor atoms to neighbors about the influence of the local environment, although it results in a high computational cost of the training due

<figure><img src="image_2.png" /><figcaption>Figure 1: 2 × 1 × 1 Pyrophyllite supercell (80 atoms), d001 represents the basal spacing. Atoms color legend: red-Oxygen, gray-Hydrogen, yellow-Silicium, purple-Aluminium.</figcaption></figure>
to the large number of parameters to fit. This method is usually more efficient with large
dataset as it helps to converge faster.

Successful MLIPs for systems similar to pyrophyllite exist in the literature. Both Kobayashi *et al.*<sup>18</sup> and Sheperd *et al.*<sup>19</sup> developed an MLIP for kaolinite clay using ACSF as descriptors. DeepMD-kit<sup>20</sup> and its smooth descriptors DeepPot-SE<sup>16</sup> have been used to train potential for various Al alloy,<sup>11</sup> mica,<sup>21</sup> and feldspar.<sup>22</sup> We also already successfully developed a MLIP for pyrophyllite clay using ACSF as descriptors within the n2p2 software.<sup>23</sup> Our representative training dataset was constructed with the well-established adaptive (or active) learning approach (allowing for great accuracy with minimal number of structures).<sup>24,25</sup>

While a dataset is usually built for a given MLIP architecture and set of descriptors, a question arises: is it possible to use the same dataset to train efficient MLIPs with different architectures and descriptors? To answer this question, we trained several MLIPs with various architectures and descriptors on an identical dataset – previously built for large-scale MD of pyrophyllite clay<sup>23</sup> – and we investigated their ability to reproduce energetic, structural, elastic and vibrationnal properties. We considered and compared an already validated MLIP using n2p2/ACSF;<sup>12,15</sup> two MLIPS built from DeePMD-kit architecture with two different types of descriptors: smooth<sup>16</sup> and attention-based descriptors;<sup>26</sup> and a message passing MLIP with learnable descriptors identical to PhysNet<sup>27</sup> but allowing to treat periodic systems.<sup>28</sup>

In addition, in another effort to assess accuracy, we investigated the effect of adding the dispersion correction *a posteriori* to the results of an MLIP fitted only to data obtained from the PBE functional. Indeed, the dispersion correction is usually already included in training datasets. Pure DFT functionals are commonly coupled with various types of corrections in order to well represent Van der Waals interactions that are crucial in numerous systems. In particular they are very important for smectite clays to well describe inter-layers interactions and it has been shown to be decisive for a proper description of silicate glasses<sup>29</sup> using MLIPs. We thus considered the widely used Grimme's D3 correction<sup>30</sup> with the BJ damping function,
based on a well-known analytical function whose contributions have been added to the n2p2
MLIP contributions. This approach's efficiency is tested by comparing with the n2p2 MLIP
fitted on PBE-D3 data.

2 Computational Details

2.1 Machine learned interaction potentials

The initial MLIP was already built from a previous study<sup>23,31</sup> using the n2p2 software<sup>32</sup> on a representative dataset of 6,974 structures obtained with adaptive learning. This dataset consists of orthorhombic and triclinic supercells computed with PBE-D3. This n2p2 MLIP is built with two hidden layer of 25 nodes each. The softplus and linear activation functions for the hidden and output layers have been used, respectively. The cut-off radius of the ACSF is set at 6 Å. For each set of radial symmetry functions, two functions were defined with the shifted radius r<sub>s</sub> kept at 0.0, two with r<sub>s</sub> shifted according to the distance of the first neighbor and two with r<sub>s</sub> shifted according to the distance of the second neighbor. The distances were taken from an equilibrium structure. Four angular functions were defined for each angle between covalently bonded atoms types.

Two DeepMD potentials<sup>20</sup>, DeepMD<sub>e2_a</sub> and DeepMD<sub>atten</sub> are trained with the two-body
embedding descriptor DeepPot-SE and with an attention-based descriptor, respectively. The
two DeepMD potentials are trained using a cut-off radius of 6 Å with a smoothing starting
at 0.5 Å for the descriptors, a 20-50-100 neural network, and a 160-160-160 embedding net.
The loss is parameterized with prefactor ρ<sub>η</sub><sup>start</sup> set at 0.02 and the prefactor ρ<sub>η</sub><sup>lim</sup> set at 1 for
the energy while ρ<sub>η</sub><sup>start</sup> is set at 1,000 and ρ<sub>η</sub><sup>lim</sup> set at 1 for the forces. The higher prefactor
on the force in the beginning is to quickly learn to fit over outliers. For the attention-based
descriptors, two attention layers are used. The training learning curves were checked to
confirm no overfitting was observed (see SI).

The PhysNet-like potential,<sup>27</sup> MPNN, is built with 5 modules containing 2 atomic residual
blocks, 3 interaction residual blocks, and 1 residual block for the output. The Gaussian basis set has 128 functions with cut-off of 5 Å. A first training is done with a ratio of 10/1 on the forces against the energy and the loss function calculated from the mean absolute error (MAE). After convergence of the loss, the ratio on the forces is increased to 20/1 and the loss calculated from RMSE to give more weight to possible outliers.

Lastly, the dataset has been recalculated with pure the PBE functional and the cor-
responding MLIP – n2p2<sub>D3</sub> – is fitted with the same parameters as the n2p2 MLIP. The
Grimme’s D3 (with the Becke-Johnson damping function correction with a cut-off of 15 Å)
is then added directly to the n2p2<sub>D3</sub> energy and forces from the dispersion/d3 pair-style
implemented in LAMMPS<sup>29</sup>.

For the training process with n2p2 and MPNN, the dataset is splited into two subsets; 90% goes to the training set and 10% goes to the test set, while for DeepMD the partitioning is 75% and 25% (standard settings were applied). The n2p2 MLIP is trained for 300 epochs, and around 1000 and 16000 epochs for MPNN and DeepMD, respectively. RMSE convergence is checked for both forces and energies and the best weights were taken along these curves (RMSE curves gathered in SI).

For all MLIPs, formatted datasets and corresponding models' files to be used directly
with LAMMPS are freely provided on Zenodo. 33

## 2.2 Simulations

The plane-wave periodic DFT calculations are made with the Vienna *ab initio* simulation package (VASP).<sup>34</sup> The PBE functional<sup>35</sup> coupled with DFT-D3 Grimme correction<sup>30</sup> is used to perform the calculations for the dataset used for the n2p2, DeepMD<sub>e2.a</sub>, DeepMD<sub>atten</sub> and MPNN MILPs, while only the PBE functional was used for the n2p2<sub>D3</sub> MLIP. It has to be specified that the D3 used for DFT calculations with VASP has been applied using the zero damping function and a cut-off of 50 Å. The sampling of the Brillouin zone is made only at the Γ-point (following the Monkhorst-Pack method). The calculations are performed with
a kinetic energy cut-off of 400 eV. These settings have already been validated in a previous work on the development of a MLIP for pyrophyllite²³ and in multiple theoretical works on smectites.³⁶,³⁷ During structure geometry optimization calculations the threshold for the force convergence is set to 0.01 eV/Å. The production for DFT MD runs are usually set to 2 ps.

The classical molecular dynamics are performed using the LAMMPS software,³⁸ using the hdnnp, DeepMD, a self developed nnmp pair style²⁸ and the hybrid dispersion/d3 – hdnnp pair styles for the n2p2, the DeepMD<sub>e2_a</sub> and DeepMD<sub>atten</sub>, the MPNN and the n2p2<sub>d3</sub> MLIPS, respectively.

### 2.2.1 Structural parameters

The structural parameters are studied after relaxation of the unit cell at 0 K and compared to experimental data. NPT simulation at 300 K, 500 K and 1,000 K are also carried out on orthorhombic and triclinic supercells to compare the averaged parameters. The radial distribution functions are obtained from a NVE trajectories on a 4 × 2 × 1 supercell after thermalization at 300 K (see SI for details).

### 2.2.2 Basal spacing study

A $d_{001}$ energy profile is obtained from MD runs on fourteen triclinic 2 × 1 × 1 supercells with a fixed $d_{001}$ ranging from 8.95 to 9.60 Å with a step size of 0.05 Å. The position of the ions on these supercells is first relaxed at the DFT level, then NVT dynamics at 300 K are performed on the relaxed cells with PBE-D3 or the MLIPs. In this case, the thermalization was performed in two runs (see SI for details). The energy profile is obtained by averaging the total energy for each production run to obtain an average energy per $d_{001}$. In addition, the averaged $d_{001}$ is computed over NPT trajectories for the unit cell, 2 × 1 × 1 and 2 × 1 × 2 supercells at 300 K, 500 K, and 1,000 K.
### 2.2.3 Exfoliation energy

The exfoliation energy, i.e. the energy required to peel off a single layer of a bulk structure, is computed by taking the difference between the unit cell energy of the bulk and that of a single layer³⁹,⁴⁰ in gas phase. It is computed as :

$$E_{exfo} = E_{bulk} - E_{layer}, \qquad (1)$$

where $E_{bulk}$ is the energy off the optimized bulk configuration and $E_{layer}$ is the energy of the optimized single layer configuration. The exfoliation energy is calculated for a c lattice vector equal to 25 Å except for MLIP$_{D3}$. Since the D3 correction has no real cut-off, the c lattice vector was increased to 100 Å to make sure there were no interactions at all.

### 2.2.4 Elastic

The elastic constant matrix $C_{ij}$ is calculated at 0 GPa on the triclinic unit cell with LAMMPS by measuring the stress tensor changes after deformation of the cell. The Voigt-Reuss-Hill approximation⁴¹ is used to calculate the bulk, shear, and elastic modulus of the pyrophyllite, as well as the Poisson's ratio.

### 2.2.5 Vibrational study

The vibrational density of state (VDoS) is evaluated from the Fourier transform of the velocity autocorrelation function obtained from NVE simulations of 10 ps at 300 K with a timestep of 1 fs on a 4 × 2 × 1 supercell.
# 3 Results and Discussion

## 3.1 Training

The RMSE values obtained during the training on the validation set are gathered in Table. 1. Both DeepMD potentials show energy RMSE values similar to the one obtained by n2p2, however, the forces RMSE values here are 2 to 3 times lower than n2p2. RMSE values close to 1 meV/atom for the energy and less than 80 meV/Å are common with DeepMD potentials for materials such as titanate, metal nanocluster, or alloy.<sup>11,20,42–44</sup> As for MPNN, the energy RMSE values obtained are 2 times larger compared to the DeepMD potentials and n2p2, however, the forces RMSE values are the lowest of our potentials. This may be due to the fact that MPNN has been trained with a loss function giving more weight to the forces compared to the energy.

## 3.2 Validation

To validate our potentials, energies and forces of each structure of an NVE *ab initio* MD of 2 ps trajectory at 300 K performed on a 2 × 1 × 2 supercell with 2 layers (configuration not included in the dataset) are compared to those computed with the MLIPs. For each step, the corresponding structure is extracted, and the energy as well as the forces on the atoms are computed using our three MLIPs. For the forces, random atoms among important atom types (surface, internal and hydroxil oxygen, hydrogen, silicium, and aluminium) are chosen and the norm of the force plotted.

The potentials all follow very closely the variation of the PBE-D3 for the energy and the

<table><caption>Table 1: RMSE values on the energy (meV/atom) and forces (meV/Å) obtained from the validation set at the end of the training of the various MLIPs.</caption><thead><tr><th></th><th>n2p2</th><th>n2p2<sub>D3</sub></th><th>DeepMD<sub>e2_a</sub></th><th>DeepMD<sub>atten</sub></th><th>MPNN</th></tr></thead><tbody><tr><th>Energy</th><td>0.9</td><td>0.9</td><td>1.6</td><td>1.7</td><td>3.0</td></tr><tr><th>Forces</th><td>183.9</td><td>182.4</td><td>66.5</td><td>86.0</td><td>40.5</td></tr></tbody></table>
forces, as can be seen in Figure. 2 and Figure. 3. The large forces fluctuations showed by the surface silicium and oxygen (Si-146 and O-55), as well as the aluminium atom (Al-120) are well-reproduced by the MLIPs. MPNN, DeepMD<sub>atten</sub>, and DeepMD<sub>e2_a</sub> show similar energy RMSE values, 0.014 meV/atom, 0.010 meV/atom, and 0.065 meV/atom, respectively, which is one order of magnitude lower compared to the MLIPs constructed with the n2p2 software (0.455 meV/atom for n2p2,<sup>23</sup> 0.533 meV/atom for n2p2<sub>D3</sub>). The same is observed for the forces RMSE, as MPNN, DeepMD<sub>atten</sub>, and DeepMD<sub>e2_a</sub> show RMSE values 2-3 times lower than the n2p2 MLIP which are around 100-150 meV/Å (Figure. 3).

n2p2 and n2p2<sub>D3</sub> MLIPs show close RMSE values for both the energy and the forces.
Adding the dispersion correction *a posteriori* to the MLIP doesn't seem to improve the
MLIP accuracy on reproducing simple MD trajectory energy and forces compared to the
MLIP fitted directly on DFT-D3.

In order to evaluate the calculation speed of each methods, a NVE MD of 2000 fs starting
from the same structure and velocities has been carried on the same CPU (See Table 2) .
The calculation times in comparison with DFT are 3 order of magnitude lower for n2p2
methods, still very efficient for DeepMD, but much less for MPNN. It has to be notices that

<figure><img src="image_3.png" /><figcaption>Figure 2: Variation of the energy on an NVE AIMD (cut-off: 400 eV) trajectory for PBE-D3, n2p2<sub>D3</sub>, MPNN, DeepMD<sub>atten</sub>, and DeepMD<sub>e2_a</sub>. The energy RMSE values are 0.533, 0.014, 0.010, and 0.065 meV/atom, respectively.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Variation of the forces on an NVE AIMD (cut-off: 400 eV) trajectory (black line), MPNN (blue line), DeepMD<sub>atten</sub> (pink line) and DeepMD<sub>e2_a</sub> (purple line). The RMSE values are in meV/Å. Si-146 is a surface silicium, O-55 is a surface oxygen atom, Al-120 is a random aluminium, O-84 is an internal oxygen, H-4 is a random hydrogen, O-66 is an hydroxil oxygen</figcaption></figure>

for the latter the mpi version is still not developped thus less intesting. However, taking into account the above validation on energy and forces and the calculation time, the MLIPs build from very well known approaches such as n2p2 and DeepMD are very beneficial compared DFT.

<table><caption>Table 2: Calculation times (in seconds) for a NVE MD of 2000 fs performed on a 64 cores computer (two CPUs AMD EPYC 7543 with 32 cores each). The initial structure and velocities were the same in all cases. The DFT calculation was performed with VASP, while NNPs MDs with LAMMPS. The MPNN module is yet not parallelized and run on 1 core.</caption><thead><tr><th></th><th>DFT</th><th>n2p2</th><th>n2p2<sub>D3</sub></th><th>DeepMD<sub>e2_a</sub></th><th>DeepMD<sub>atten</sub></th><th>MPNN</th></tr></thead><tbody><tr><td>Calc. time</td><td>6168.0</td><td>3.2</td><td>6.3</td><td>16.6</td><td>36.0</td><td>491.3</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 4: Energy profile of the basal spacing ($d_{001}$) calculated with PBE-D3, n2p2, n2p2$_{D3}$, DeepMD$_{e2\_a}$, DeepMD$_{atten}$, and MPNN.</figcaption></figure>

## 3.3 Energy profile of the $d_{001}$

Following our validation protocol, the profile of the energy evolution according to the $d_{001}$ is investigated. Figure 4 shows the energy profile of the basal spacing obtained from averaged values computed on NVT runs at 300 K. The PBE-D3 curve shows a minimum at 9.30 Å, n2p2 at 9.40 Å. The DeepMD$_{e2\_a}$ and DeepMD$_{atten}$ curves minimum are at 9.35 Å and 9.40 Å, respectively, showing a flatter profile. These MLIPs all describe quite well the repulsive regime, and follow the PBE-D3 curve closely, however, differences start to show in the attractive regime. The MPNN curve starts to slightly diverge and overestimate the energy after 9.30 Å seemingly giving more weight to attraction, which might be due to the lower cut-off used for its descriptors. The DeepMD potentials and MPNN curves are rather close until the $d_{001}$ reaches 9.30 Å, at which point they start to very slightly diverge; the difference between the DeepMD$_{atten}$ and DeepMD$_{e2\_a}$ curves is negligible. On the contrary, n2p2 has a tendency to underestimate the energy in this region, being more accurate in the repulsive regime. n2p2$_{D3}$ surprisingly fails to reproduce the proper PBE-D3 behavior in the attractive regime. This could be due to a cut-off issue as the MLIP interactions are
subjected to the descriptor's cut-off while the D3 correction doesn't, thus leading to some discontinuities.

Simulations in the NPT ensemble at 300 K, 500 K and 1,000 K are run on various pyrophyllite supercells and the averaged $d_{001}$ values are compared with PBE-D3 results (see Table. 3). For the DFT calculations on the unit cell, one can observe a very weak decrease in $d_{001}$ from 300 to 500 or 1000 K of about 0.032 and 0.018 Å. However, these results should be considered in light of the fact that they relate to a very small cell containing only 40 atoms, which means they reflect thermal expansion less precisely. All MLIPs show results rather close to PBE-D3, although they all overestimate the $d_{001}$ value, and the MLIPs accuracy in predicting the $d_{001}$ slightly decreases when increasing the temperature. Moreover, the results are more accurate for the larger cells (2 × 1 × 1 and 2 × 1 × 2) than for the unit cell. This is especially true for n2p2, which seems better at extrapolating to larger cells than to smaller ones. At 300 K and 500 K, the best results are obtained with n2p2 for the three cells, while

<table><caption>Table 3: Averaged $d_{001}$ values (in Å) obtained from NPT runs for different supercells sizes and different temperatures.</caption><thead><tr><th></th><th>unit cell</th><th>2 × 1 × 1</th><th>2 × 1 × 2</th></tr></thead><tbody><tr><td colspan="4"><strong>300 K</strong></td></tr><tr><td>PBE-D3<sup>a</sup></td><td>9.287 ± 0.070</td><td>9.257 ± 0.060</td><td>9.298 ± 0.020</td></tr><tr><td>n2p2</td><td>9.431 ± 0.030</td><td>9.427 ± 0.003</td><td>9.422 ± 0.003</td></tr><tr><td>n2p2<sub>D3</sub></td><td>9.453 ± 0.028</td><td>9.434 ± 0.033</td><td>9.437 ± 0.022</td></tr><tr><td>DeepMD<sub>e2.a</sub></td><td>9.538 ± 0.041</td><td>9.445 ± 0.011</td><td>9.451 ± 0.010</td></tr><tr><td>DeepMD<sub>atten</sub></td><td>9.576 ± 0.042</td><td>9.506 ± 0.029</td><td>9.489 ± 0.011</td></tr><tr><td colspan="4"><strong>500 K</strong></td></tr><tr><td>PBE-D3<sup>a</sup></td><td>9.255 ± 0.100</td><td>9.342 ± 0.100</td><td>9.333 ± 0.040</td></tr><tr><td>n2p2</td><td>9.504 ± 0.030</td><td>9.438 ± 0.006</td><td>9.433 ± 0.005</td></tr><tr><td>n2p2<sub>D3</sub></td><td>9.517 ± 0.116</td><td>9.435 ± 0.044</td><td>9.442 ± 0.009</td></tr><tr><td>DeepMD<sub>e2.a</sub></td><td>9.675 ± 0.027</td><td>9.508 ± 0.032</td><td>9.573 ± 0.015</td></tr><tr><td>DeepMD<sub>atten</sub></td><td>9.712 ± 0.036</td><td>9.564 ± 0.025</td><td>9.575 ± 0.025</td></tr><tr><td colspan="4"><strong>1,000 K</strong></td></tr><tr><td>PBE-D3<sup>a</sup></td><td>9.269 ± 0.100</td><td>9.527 ± 0.100</td><td>9.430 ± 0.080</td></tr><tr><td>n2p2</td><td>9.962 ± 0.300</td><td>9.803 ± 0.200</td><td>9.686 ± 0.070</td></tr><tr><td>n2p2<sub>D3</sub></td><td>9.636 ± 0.224</td><td>9.483 ± 0.050</td><td>9.556 ± 0.668</td></tr><tr><td>DeepMD<sub>e2.a</sub></td><td>9.782 ± 0.101</td><td>9.723 ± 0.031</td><td>9.749 ± 0.026</td></tr><tr><td>DeepMD<sub>atten</sub></td><td>9.857 ± 0.074</td><td>9.761 ± 0.041</td><td>9.742 ± 0.016</td></tr></tbody></table>

<footnote>a energy cutoff of 600 eV used for DFT calculations.</footnote>
at 1,000 K the two DeepMD potentials show more accurate results – especially DeepMD<sub>e2_a</sub>. It thus appears that the n2p2 potential is better at describing equilibrium structures, while the DeepMD potentials extrapolate better as they are being less affected by the increased temperature. The better extrapolation capacity of the DeepMD model might be explained by two points: i) the slightly more complex architecture (with more neurons and hidden layers) of the DeepMD NN compared to our n2p2 NN; and ii) DeepMD descriptors include information about the local atomic environment in a more comprehensive and systematic way than the ACSF. n2p2<sub>D3</sub> gives results close to PBE-D3 *d*<sub>001</sub> even at 1,000 K with the larger overestimation being only of 0.367 Å (unit cell at 1,000 K). n2p2<sub>D3</sub> also gives stable results when increasing cell size and temperature. The difference between the PBE-D3 and n2p2 MLIP results increases from 0.1 Å to 0.3 Å for the 2 × 1 × 2 supercell when the temperature increases from 500 K to 1,000 K, while the difference with the n2p2<sub>D3</sub> results stay at 0.1 Å. This shows the stabilizing effect of adding the D3 correction *a* *posteriori* to the MLIP (which is not sensitive to the MLIP's cutoff effect). However, this is not always the case, as the *d*<sub>001</sub> energy profile shows the n2p2<sub>D3</sub>'s poor accuracy for larger values than 9.3 Å.

## 3.4 Structural analysis

The unit cell and ionic positions are relaxed at 0 K and the lattice parameters are listed in Table. 4. The results for all MLIPs are very close to the PBE-D3 reference, with an error

<table><caption>Table 4: Lattice parameters after relaxation of the unit cell and ionic positions (errors to PBE-D3 value between parentheses). Distances: *a*, *b*, *c* and *d*<sub>001</sub> given in Å, angles in degrees.</caption><thead><tr><th>0 K relaxation</th><th>a</th><th>b</th><th>d<sub>001</sub></th><th>α</th><th>β</th><th>γ</th><th>volume</th></tr></thead><tbody><tr><td>PBE-D3<sup>a</sup></td><td>5.134</td><td>8.921</td><td>9.300</td><td>90.7</td><td>100.0</td><td>89.2</td><td>432.3</td></tr><tr><td>n2p2</td><td>5.202 (+1.3%)</td><td>9.059 (+1.5%)</td><td>9.319 (+0.2%)</td><td>91.7 (+1.1%)</td><td>100.5 (+0.5%)</td><td>90.1 (+1.0%)</td><td>439.2 (+1.6%)</td></tr><tr><td>n2p2<sub>D3</sub></td><td>5.175 (+0.8%)</td><td>9.010 (+1.0%)</td><td>9.318 (+0.2%)</td><td>91.9 (+1.3%)</td><td>100.7 (+0.7%)</td><td>89.9 (+0.8%)</td><td>434.4 (+0.5%)</td></tr><tr><td>DeepMD<sub>e2_a</sub></td><td>5.193 (+1.1%)</td><td>9.004 (+0.9%)</td><td>9.352 (+0.5%)</td><td>90.6 (-0.1%)</td><td>100.8 (+0.8%)</td><td>90.0 (+0.9%)</td><td>437.3 (+1.2%)</td></tr><tr><td>DeepMD<sub>atten</sub></td><td>5.199 (+1.3%)</td><td>9.018 (+1.1%)</td><td>9.410 (+1.2%)</td><td>90.7 (0.0%)</td><td>100.9 (+0.9%)</td><td>89.9 (+0.8%)</td><td>441.3 (+2.1%)</td></tr></tbody></table>

<footnote><sup>a</sup> energy cutoff of 400 eV used for DFT calculations.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5: Data of $g(r)$ in Å obtained from an NVE run at 300 K on a 4 × 2 supercell. a) Shift of the first neighbor peak for various atom pairs. b) $g(r)$ computed for the Si–Si pair. c) $R^2$ score of the $g(r)$ plot with respect to PBE-D3 reference. d) $g(r)$ computed for the H–O pair.</figcaption></figure>

of less than 3%. DeepMD<sub>atten</sub> shows higher error for the $d_{001}$ value, 1.2% against 0.9% and 0.5% for DeepMD<sub>e2_a</sub> and n2p2 respectively, and the volume estimation with an error of 2.1% against 1.2% and 1.6% for DeepMD<sub>e2_a</sub> and n2p2 respectively. n2p2<sub>D3</sub> is overall the most accurate with respect to PBE-D3, except for the angle $\alpha$ where it shows a slightly higher error of +1.3%. Moreover, its error on the volume is the smallest, only 0.5%, showing that the error on the angle $\alpha$ should mostly be due to more important layers sliding during MLIP<sub>D3</sub> simulations.

The radial distribution functions $g(r)$ are analysed along an NVE production run. The shift of the first peak compared to PBE-D3 and the $R^2$ score of the function as well as two
example curves of the Si-Si pairs and H-O pairs are gathered in Figure. 5. Globally, as seen in Figure. 5 a), the first peak shift is quite small, inferior to 0.05 Å for all MLIPs, thus the four MLIPs reproduce extremely well the position of the first neighbor peak of the radial distribution functions. n2p2, n2p2<sub>D3</sub> and DeepMD<sub>e2_a</sub> have similar shift for all pairs. This is not surprising as both ACSF and DeepMD-kit's DeepPot-SE descriptors contain directly and explicitly radial information of the local environment. One can observe for the Si-Si pair (Figure. 5b) that the second peak is slightly over estimated for DeepMD and n2p2 MLIPs, of about 0.2-3 Å. Here even the more comprehensive description of DeepMD descriptors, the embedding attention-based ones, slightly fail to accurately describe the behaviour of the second neighbors shell. In comparison, the message passing descriptors based MLIP seem to perform better. The MLIPs show good R<sup>2</sup> score for most atom pairs with especially good score for O-Al, O-Si, and O-O pairs. For O-H the bad R<sup>2</sup> score for MPNN is only due to a higher, more thin first peak of g(r).

## 3.5 Exfoliation

The exfoliation energies are gathered in Table. 5. Both DeepMD<sub>e2_a</sub> and DeepMD<sub>atten</sub> are close to the PBE-D3 exfoliation energy with differences of 0.07 eV and 0.096 eV, respectively. Although, the exfoliation energy predicted by the MPNN is the closest to the PBE-D3 result with only a small difference of 0.023 eV. n2p2 shows the larger differences of 0.113 eV to the PBE-D3 exfoliation energy compared to the other MLIPs, this is consistent with the fact that the energy RMSE values of the validation MD trajectory (section 3.2) are one order of magnitude lower for these three potentials compared to n2p2. Adding the dispersion correction a posteriori gives better results than including it in the training set.

<table><caption>Table 5: Pyrophyllite exfoliation energy E<sub>exfo</sub> in eV, values between parenthesis are in kJ.mol<sup>-1</sup> for comparison purpose.</caption><thead><tr><th></th><th>PBE-D3</th><th>n2p2</th><th>n2p2<sub>D3</sub></th><th>MPNN</th><th>DeepMD<sub>e2_a</sub></th><th>DeepMD<sub>atten</sub></th></tr></thead><tbody><tr><th>E<sub>exfo</sub></th><td>-0.621 (-60)</td><td>-0.508 (-49)</td><td>-0.574 (-56)</td><td>-0.597 (-58)</td><td>-0.551 (-53)</td><td>-0.525 (-51)</td></tr></tbody></table>
The difference with PBE-D3 decreases from 0.113 eV for n2p2 to 0.047 eV for n2p2<sub>D3</sub>. Indeed, as the pyrophyllite layer are held through van der Waals, an accurate treatment of the dispersion correction is crucial to correctly model the exfoliation.

## 3.6 Elastic modulii.

Elastic modulii are gathered in Table. 6. Calculating these values at the DFT level requires significant computer resources. Therefore, MLIP values are compared to DFT values computed in previous work by Qin et al.,<sup>46</sup> which used the slightly different D2 dispersion correction. All MLIPs overestimate the Young's modulus E, with predicted values close to double the PBE-D2 value. Although, n2p2 is the closest with 37 GPa of difference while DeepMD<sub>e2_a</sub> is the farthest with a difference of 47 GPa with respect to the PBE-D2 result. They also all slightly overestimate the shear modulus G, with predicted values close to double the PBE-D2 value. n2p2 being the closest to the PBE-D2 result with a difference of 17 GPa while DeepMD<sub>e2_a</sub> is the least accurate one with a difference of 28 GPa. However, the bulk modulus B values of the MLIPs are quite close to the PBE-D2 result and the experimental range of values. DeepMD<sub>e2_a</sub> is the more accurate MLIPs with a difference of only 1.2 GPa, DeepMD<sub>atten</sub> is also quite close with a difference of 1.6 GPa. Regarding the Poisson's ratio, n2p2 yields a more accurate result (to the PBE-D2 results and experimental range of values) than the two DeepMD potentials. Overall, the differences in accuracy between the three MLIPs are minimal.

<table><caption>Table 6: Bulk modulus (B), shear modulus (G), Young's modulus (E) and Poisson's ratio (ν) in GPa obtained from calculations at 0 GPa with the developed MLIPs and compared to different studies.</caption><thead><tr><th></th><th>Experiments<sup>45</sup></th><th>PBE-D2<sup>a</sup></th><th>n2p2</th><th>n2p2<sub>D3</sub></th><th>DeepMD<sub>e2_a</sub></th><th>DeepMD<sub>atten</sub></th></tr></thead><tbody><tr><td>B</td><td>14.80–37.50</td><td>44.64</td><td>39.57</td><td>41.05</td><td>40.80</td><td>37.96</td></tr><tr><td>G</td><td>18.20–21.30</td><td>27.47</td><td>44.67</td><td>49.40</td><td>55.88</td><td>51.63</td></tr><tr><td>E</td><td>38.60–53.60</td><td>68.39</td><td>104.02</td><td>105.78</td><td>115.09</td><td>106.58</td></tr><tr><td>ν</td><td>0.06–0.26</td><td>0.25</td><td>0.16</td><td>0.07</td><td>0.03</td><td>0.03</td></tr></tbody></table>

<footnote><sup>a</sup> The PBE-D2 values are taken from ref: Qin *et al.*<sup>46</sup></footnote>
The differences between the results produced by n2p2 and n2p2<sub>D3</sub> are minimal as well. Their results differ by only around 2 GPa for the bulk modulus B and the Young's modulus E, and by at most 5 GPa for the shear modulus G.

## 3.7 Vibrational analysis

All the MLIPs show rather similar spectra, with one peak for the OH bending around 950 cm⁻¹, peaks for the Si-O vibrations in the 1000–1200 cm⁻¹ region, and O-H stretching peaks around 3800 cm⁻¹. However, the PBE-D3 spectrum in the 1000–1150 cm⁻¹ region (associated with Si-O vibrations⁴⁷) differs from the MLIPs. PBE-D3 only seems to show a low intensity peak at 1020 cm⁻¹ while the MLIPs clearly show Si-O vibration peaks around 1010–1030 cm⁻¹ and 1080–1150 cm⁻¹.

The OH bending peaks showed by DeepMD<sub>e2_a</sub>, DeepMD<sub>atten</sub>, and MPNN are less intense than the peaks from the n2p2 and n2p2<sub>D3</sub> and even much less intense than the PBE-D3 peak. Although, the MLIPs OH bending peaks positions are quite accurate to PBE-D3, located at around 920 cm⁻¹ for all MLIPs with a shift of only 30 cm⁻¹ compared to the PBE-D3

<figure><img src="image_7.png" /><figcaption>Figure 6: Vibrational density of state of pyrophyllite computed from 10 ps MDs with PBE-D3, n2p2, n2p2<sub>D3</sub>, DeepMD<sub>e2_a</sub>, DeepMD<sub>atten</sub>, MPNN, experimental data taken from.⁴⁷</figcaption></figure>
spectrum. The PBE-D3 OH stretching double peak is located at 3795-3805 cm⁻¹. The MLIPs reproduce well this peak with minimal shift; DeepMD<sub>e2_a</sub> and DeepMD<sub>atten</sub> show a single peak placed at 3800 cm⁻¹; MPNN also shows a single peak at 3810 cm⁻¹. n2p2 and n2p2<sub>D3</sub> show two peaks at 3775 cm⁻¹ and 3810 cm⁻¹, however, n2p2<sub>D3</sub> double peak is very weak.

The DeepMD potentials (DeepMD<sub>e2_a</sub> and DeepMD<sub>atten</sub>) and MPNN show extremely similar spectra with peaks at the same positions with only slight differences in the intensities of the peaks. The n2p2 and n2p2<sub>D3</sub> spectra are also very similar to each other. Although, it can be noted that these two spectra show peak forms, positions and intensities that are slightly different from the other three potentials (DeepMD<sub>e2_a</sub>, DeepMD<sub>atten</sub> and MPNN). Overall, all MLIPs show OH bending and stretching peaks close to the ones in the PBE-D3 spectrum.

# 4 Conclusion

Different MLIPs have been trained on an identical representative dataset obtained by an adaptive learning procedure. Two potentials have been trained with the DeepMD-kit package, one with a two-body embedding descriptor (DeepMD<sub>e2_a</sub>), and the other with an attention-based descriptor (DeepMD<sub>atten</sub>). A third potential was created using a message-passing graph neural network based on the PhysNet model (MPNN). These three potentials have been compared to the performances of the previous n2p2, build with ACSF as descriptor, on their accuracy to the PBE-D3 reference calculations. The portability of the dataset – the ability of the dataset to yield accurate potentials when trained with different models – as well as the possible influence of the different descriptors on the accuracy of the potentials was investigated. To better describe the dispersion forces, a MLIP was trained on pure PBE results with Grimme's D3 correction added *a posteriori* to avoid cut-off. All the MLIPs were able to reproduce PBE-D3 energies and forces very accurately, which is the most important
to obtain accurate MD trajectories.

The MLIP with the dispersion correction added *a posteriori* perform more accurately
than its counterpart trained on PBE-D3 (n2p2). Especially for structural parameters where
n2p2<sub>D3</sub> shows more stable and accurate results over various temperature and cell sizes.
n2p2<sub>D3</sub> shows also a more accurate exfoliation energy.

The four MLIPs trained on the same dataset (n2p2, DeepMD<sub>e2_a</sub>, DeepMD<sub>atten</sub>, MPNN)
were very close in accuracy for the various parameters investigated: structural parameters,
<i>d</i><sub>001</sub> energy profile, exfoliation energy, elastic properties, and vibrational spectra. One can
see that if the dataset is already representative, changing the descriptors doesn't have a lot
of influence on the accuracy of the obtained potentials.

**Supporting Information Available**

Supporting information contains details of the MD simulations at the DFT level as well as
performed with the MLIPs. In addition, the energy and forces training convergence curves
on the are included for each MLIP.

Acknowledgement

C.S and P.M thank the CNRS support for granting this project through the program
INSU-EC2CO/HYBIIGE (project RENOIR) and the program 80|PRIME-2022 (project Pi-
casso). This work was performed using HPC resources from GENCI-IDRIS (Grant A17-
A0170807662). We gratefully acknowledge support from the CNRS/IN2P3 Computing Cen-
ter (Lyon - France) for providing computing and data-processing resources needed for this
work.
References

(1) Bergaya, F.; Lagaly, G. *Developments in Clay Science*; Elsevier, 2006; Vol. 1; Chapter 1 General Introduction: Clays, Clay Minerals, and Clay Science, pp 1–18, DOI: doi:10.1016/S1572-4352(05)01001-9.

(2) Vaccari, A. Clays and catalysis: a promising future. *Applied Clay Science* **1999**, *14*, 161–198, DOI: doi:10.1016/S0169-1317(98)00058-1.

(3) Mccabe, R.; Adams, J. *Developments in Clay Science*; Elsevier, 2013; Vol. 5; Chapter 4.3 Clay Minerals as Catalysts, pp 491–538, DOI: doi:10.1016/B978-0-08-098259-5.00019-6.

(4) Ochirkhuyag, A.; Temuujin, J. The Catalytic Potential of Modified Clays: A Review. *Minerals* **2024**, *14*, 629, DOI: doi:10.3390/min14060629.

(5) Bergaya, F.; Theng, B.; Lagaly, G. In *Handbook of Clay Science*; Bergaya, F., Theng, B. K., Lagaly, G., Eds.; Developments in Clay Science; Elsevier, 2006; Vol. 1; Chapter 10 Clays in Industry, p 499, DOI: doi:https://doi.org/10.1016/S1572-4352(05)01015-9.

(6) Bergaya, F.; Theng, B.; Lagaly, G. In *Handbook of Clay Science*; Bergaya, F., Theng, B. K., Lagaly, G., Eds.; Developments in Clay Science; Elsevier, 2006; Vol. 1; Chapter 11 Clays, Environment and Health, p 623, DOI: doi:https://doi.org/10.1016/S1572-4352(05)01019-6.

(7) Thiebault, T. Raw and modified clays and clay minerals for the removal of pharmaceutical products from aqueous solutions: State of the art and future perspectives. *Critical Reviews in Environmental Science and Technology* **2020**, *50*, 1451–1514, DOI: doi:10.1080/10643389.2019.1663065.
(8) Zhang, J.; Zhou, Y.; Jiang, M.; Li, J.; Sheng, J. Removal of methylene blue from aqueous solution by adsorption on pyrophyllite. *Journal of Molecular Liquids* **2015**, *209*, 267–271, DOI: doi:10.1016/j.molliq.2015.05.056.

(9) Hasanbegović, E.; Huremović, J.; Žero, S. Adsorption capacity of nitrate from artificial fertilizers and soil on pyrophyllite. *International Journal of Environmental Science and Technology* **2021**, *18*, 3731–3738, DOI: doi:10.1007/s13762-021-03135-2.

(10) Qin, X.; Zhao, J.; Wang, J.; He, M. Atomic Structure, Electronic and Mechanical Properties of Pyrophyllite under Pressure: A First-Principles Study. *Minerals* **2020**, *10*, 778, DOI: doi:10.3390/min10090778.

(11) Wen, T.; Zhang, L.; Wang, H.; E, W.; Srolovitz, D. J. Deep potentials for materials science. *Materials Futures* **2022**, *1*, 022601, DOI: doi:10.1088/2752-5724/ac681d.

(12) Behler, J. Atom-centered symmetry functions for constructing high-dimensional neural network potentials. *The Journal of Chemical Physics* **2011**, *134*, 074106, DOI: doi:10.1063/1.3553717.

(13) Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. *Phys. Rev. B* **2013**, *87*, 184115, DOI: doi:10.1103/PhysRevB.87.184115.

(14) Drautz, R. Atomic cluster expansion for accurate and transferable interatomic potentials. *Phys. Rev. B* **2019**, *99*, 014104, DOI: doi:10.1103/PhysRevB.99.014104.

(15) Behler, J.; Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. *Physical Review Letters* **2007**, *98*, 146401, DOI: doi:10.1103/PhysRevLett.98.146401.

(16) Zhang, L.; Han, J.; Wang, H.; Saidi, W. A.; Car, R.; Weinan, E. End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems.
*Proceedings of the 32nd International Conference on Neural Information Processing Systems* **2018**, 4441–4451, DOI: doi:10.48550/ARXIV.1805.09003.

(17) Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; Dahl, G. E. Neural message passing for Quantum chemistry. *Proceedings of the 34th International Conference on Machine Learning - Volume 70* **2017**, 1263–1272, DOI: doi:10.48550/ARXIV.1704.01212.

(18) Kobayashi, K.; Yamaguchi, A.; Okumura, M. Machine learning potentials of kaolinite based on the potential energy surfaces of GGA and meta-GGA density functional theory. *Applied Clay Science* **2022**, 228, 106596, DOI: doi:10.1016/j.clay.2022.106596.

(19) Shepherd, S.; Tribello, G. A.; Wilkins, D. M. A Fully Quantum-Mechanical Treatment for Kaolinite. *The Journal of Chemical Physics* **2023**, 158, 204704, DOI: doi:10.1063/5.0152361.

(20) Zeng, J.; Zhang, D.; Lu, D.; Mo, P.; Li, Z.; Chen, Y.; Rynik, M.; Huang, L.; Li, Z.; Shi, S.; Wang, Y.; Ye, H.; Tuo, P.; Yang, J.; Ding, Y.; Li, Y.; Tisi, D.; Zeng, Q.; Bao, H.; Xia, Y.; Huang, J.; Muraoka, K.; Wang, Y.; Chang, J.; Yuan, F.; Bore, S. L.; Cai, C.; Lin, Y.; Wang, B.; Xu, J.; Zhu, J.-X.; Luo, C.; Zhang, Y.; Goodall, R. E. A.; Liang, W.; Singh, A. K.; Yao, S.; Zhang, J.; Wentzcovitch, R.; Han, J.; Liu, J.; Jia, W.; York, D. M.; E, W.; Car, R.; Zhang, L.; Wang, H. DeePMD-kit v2: A software package for deep potential models. *The Journal of Chemical Physics* **2023**, 159, 054801, DOI: doi:10.1063/5.0155600.

(21) Raman, A. S.; Selloni, A. An *ab-initio* deep neural network potential for accurate large-scale simulations of the muscovite mica-water interface. *Molecular Physics* **2025**, 123, e2365430, DOI: doi:10.1080/00268976.2024.2365430.

(22) Piaggi, P. M.; Selloni, A.; Panagiotopoulos, A. Z.; Car, R.; Debenedetti, P. G. A first-principles machine-learning force field for heterogeneous ice nucleation on microcline feldspar. *Faraday Discussions* **2024**, 249, 98–113, DOI: doi:10.1039/D3FD00100H.
(23) Sanz, C.; Allouche, A.-R.; Bousige, C.; Mignon, P. Neural Network Atomistic Potential for Pyrophyllite Clay Simulations. *The Journal of Physical Chemistry A* **2025**, *129*, 3567–3577, DOI: doi:10.1021/acs.jpca.5c00406.

(24) Mignon, P.; Allouche, A.-R.; Innis, N. R.; Bousige, C. Neural Network Approach for a Rapid Prediction of Metal-Supported Borophene Properties. *Journal of the American Chemical Society* **2023**, *145*, 27857–27866, DOI: doi:10.1021/jacs.3c11549.

(25) Bousige, C.; Delenda, A.-A.; Allouche, A.-R.; Mignon, P. A Portable Data Set for Borophene Growth Modeling with Reactive Neural Network Potentials. *J. Phys. Chem. C* **2025**, *129*, 18760–18771, DOI: doi:10.1021/acs.jpcc.5c04912.

(26) Zhang, D.; Bi, H.; Dai, F.-Z.; Jiang, W.; Zhang, L.; Wang, H. Pretraining of attention-based deep learning potential model for molecular simulation. *npj Computational Materials* **2024**, *94*, DOI: doi:10.1038/s41524-024-01278-7.

(27) Unke, O. T.; Meuwly, M. PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges. *Journal of Chemical Theory and Computation* **2019**, *15*, 3678–3693, DOI: doi:10.1021/acs.jctc.9b00181.

(28) Allouche, A. Allouchear/NNMP-Pot. 2024 (accessed 2025-03-04); https://github.com/allouchear/NNMP-Pot.

(29) Pedone, A.; Bertani, M.; Benassi, M. Revisiting Machine Learning Potentials for Silicate Glasses: The Missing Role of Dispersion Interactions. *Journal of Chemical Theory and Computation* **2025**, *21*, 4769–4778, DOI: doi:10.1021/acs.jctc.5c00218, PMID: 40272457.

(30) 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. *The Journal of Chemical Physics* **2010**, *132*, 154104, DOI: doi:10.1063/1.3382344.
(31) Mignon, P. Datasets and Models of Pyropyllite Machine Learned Interaction Potentials, for n2p2, DeepMD_e2, DeepMD attention, PhysNet Potential types. 2025; https://doi.org/10.5281/zenodo.17395141.

(32) Singraber, A.; Bircher, M.; Reeve, S.; Swenson, D. W. H.; Lauret, J.; David, P. Comp-
PhysVienna/N2p2: Version 2.1.4. Zenodo 2021, DOI: doi:10.5281/zenodo.4750573.

(33) Mignon, P. Datasets of DFT/PBE structures used by neural network atomistic potentials for Pyrophyllite smectite Clay. 2025; https://doi.org/10.5281/zenodo.14608336.

(34) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Physical Review B* **1996**, *54*, 11169–11186, DOI: doi:10.1103/PhysRevB.54.11169.

(35) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. *Physical Review Letters* **1996**, *77*, 3865–3868, DOI: doi:10.1103/PhysRevLett.77.3865.

(36) Mignon, P.; Ugliengo, P.; Sodupe, M.; Hernandez, E. R. Ab initio molecular dynamics study of the hydration of Li+, Na+ and K+ in a montmorillonite model. Influence of isomorphic substitution. *Phys. Chem. Chem. Phys.* **2010**, *12*, 688–697, DOI: doi:10.1039/B915689E.

(37) Mignon, P.; Sodupe, M. Structural Behaviors of Cytosine into the Hydrated Interlayer of Na+-Montmorillonite Clay. An ab Initio Molecular Dynamics Study. *The Journal of Physical Chemistry C* **2013**, *117*, 26179–26189, DOI: doi:10.1021/jp4103383.

(38) 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.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a
flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Computer Physics Communications* **2022**, *271*, 108171, DOI: doi:10.1016/j.cpc.2021.108171.

(39) Jung, J. H.; Park, C.-H.; Ihm, J. A Rigorous Method of Calculating Exfoliation Energies from First Principles. *Nano Letters* **2018**, *18*, 2759-2765, DOI: doi:10.1021/acs.nanolett.7b04201.

(40) Cutini, M.; Maschio, L.; Ugliengo, P. Exfoliation Energy of Layered Materials by DFT-D: Beware of Dispersion! *Journal of Chemical Theory and Computation* **2020**, *16*, 5244-5252, DOI: doi:10.1021/acs.jctc.0c00149.

(41) Hill, R. The Elastic Behaviour of a Crystalline Aggregate. *Proceedings of the Physical Society. Section A* **1952**, *65*, 349-354, DOI: doi:10.1088/0370-1298/65/5/307.

(42) Gao, Y.; Wang, X.; Yuan, H.; Xu, H. Thermal conductivity of the layered titanate K0.8Li0.27Ti1.73O4 explored by a deep learning interatomic potential. *The Journal of Chemical Physics* **2025**, *162*, 124702, DOI: doi:10.1063/5.0255515.

(43) Tiwari, V.; Karmakar, T. Atomistic Details of Nanocluster Formation from Machine-Learned-Potential-Based Simulations. *Nano Letters* **2025**, *25*, 5940-5946, DOI: doi:10.1021/acs.nanolett.5c01340.

(44) Jiang, W.; Zhang, Y.; Zhang, L.; Wang, H. Accurate Deep Potential model for the Al-Cu-Mg alloy in the full concentration space*. *Chinese Physics B* **2021**, *30*, 050706, DOI: doi:10.1088/1674-1056/abf134.

(45) Aleksandrov, K. S.; Prodaivoda, G. T. Elastic properties of minerals. *Crystallography Reports* **1993**, *38*.

(46) Qin, X.; Zhao, J.; Wei, R.; Xu, X.; He, M. Mechanical properties of pyrophyllite under
the coupling of high temperature-pressure: A first-principle study. *Applied Clay Science* **2022**, *228*, 106613, DOI: doi:10.1016/j.clay.2022.106613.

(47) Farmer, V.; Russell, J. The infra-red spectra of layer silicates. *Spectrochimica Acta* **1964**, *20*, 1149-1173, DOI: doi:https://doi.org/10.1016/0371-1951(64)80165-X.
<figure><img src="image_8.png" /><figcaption>Figure 7: TOC.</figcaption></figure>