# Structure and Properties of Water in a New Model of the 10-Å Phase: Classical and Ab Initio Atomistic Computational Modeling

Evgeny V Tararushkin, Grigory S Smirnov, Andrey G. Kalinichev

► To cite this version:

Evgeny V Tararushkin, Grigory S Smirnov, Andrey G. Kalinichev. Structure and Properties of Water in a New Model of the 10-Å Phase: Classical and Ab Initio Atomistic Computational Modeling. Minerals, 2023, 13 (8), pp.1018. 10.3390/min13081018. hal-04197423

HAL Id: hal-04197423

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

Submitted on 15 May 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Structure and Properties of Water in a New Model of the 10-Å Phase: Classical and Ab Initio Atomistic Computational Modeling

Evgeny V. Tararushkin<fnref n="1" />, Grigory S. Smirnov<fnref n="1" /> and Andrey G. Kalinichev<fnref n="2" />,*

**Abstract:** The 10-Å phase is an important member of the family of dense hydrous magnesium silicates (DHMSs) that play a major role in the water budget in the Earth's upper mantle. Its nominal composition is usually written as $Mg_3Si_4O_{10}(OH)_2 \cdot xH_2O$, and its structure is often described as layers of talc with some amount of water present in the interlayer space. However, its actual structure and composition and the detailed mechanisms of retaining $H_2O$ molecules within the mineral are not yet sufficiently known. In particular, more recent spectroscopic and diffraction data indicate the presence of Si vacancies in the tetrahedral silicate sheets of the 10-Å phase leading to the formation of $Q^2$-type Si sites terminated by silanol groups. These silanols are, in turn, hydrogen bonded to interlayer $H_2O$ molecules. Here, we use classical and ab initio molecular dynamics (MD) simulations to compare the structures and properties of ideal and defect models of the 10-Å phase under ambient conditions. For classical MD simulations, the most recent modification of the ClayFF force field is used, which can accurately account for the bending of Mg–O–H and Si–O–H angles in the mineral layers, including the structural defects. The crystal lattice parameters, elastic constants, structure, and dynamics of the interlayer hydrogen bonding network for the model 10-Å phase are calculated and compared with available experimental data. The results demonstrate that the inclusion of Si vacancies leads to better agreement with crystallographic data, elastic constants, and bulk and shear moduli compared to a simpler model based on the idealized talc structure. The results also clearly illustrate the importance of the explicit inclusion of Mg–O–H and Si–O–H angular bending terms for accurate modeling of the 10-Å phase. In particular, the properly constrained orientation of the silanol groups promotes the formation of strong hydrogen bonds with the interlayer $H_2O$ molecules.

**Keywords:** water in minerals; talc; 10-Å phase; molecular dynamics simulations; ab initio calculations; crystal structure; elastic constants; vibrational spectra; hydrogen bonding

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

### check for updates

**Citation:** Tararushkin, E.V.; Smirnov, G.S.; Kalinichev, A.G. Structure and Properties of Water in a New Model of the 10-Å Phase: Classical and Ab Initio Atomistic Computational Modeling. *Minerals* **2023**, *13*, 1018. https://doi.org/10.3390/min13081018

Academic Editors: Feiwu Zhang, Jin Liu and Joshua Muir

Received: 13 June 2023

Revised: 19 July 2023

Accepted: 27 July 2023

Published: 30 July 2023

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

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The so-called 10-Å phase (TAP) belongs to the family of dense hydrous magnesium silicates (DHMS). It plays a key role in transporting and storing water in the Earth's mantle at subduction zones [1–5]. Therefore, TAP structural, thermodynamic, and mechanical properties are important, especially at high pressures and temperatures [6–9]. TAP is assumed to have the composition of $Mg_3Si_4O_{10}(OH)_2 \cdot xH_2O$ and consists of talc-like T–O–T layers where one sheet (O) of octahedrally coordinated Mg atoms is sandwiched from both sides by two sheets (T) of tetrahedrally coordinated Si atoms. In contrast to hydrophobic talc, TAP is assumed to contain some amount of $H_2O$ molecules in the interlayer between its T–O–T layers. Water contents of $x = 2/3$ [7], $x = 1.0$ [10], and $x = 2.0$ [6] are most commonly assumed (Figure 1). The $H_2O$ molecules occupy the six-membered siloxane rings of the tetrahedral sheets. According to experimental results, TAP has different types of symmetry [9,11]. Comodi et al. [11] showed that the TAP structure is very similar to that

<footnote n="1">¹ International Laboratory for Supercomputer Atomistic Modelling and Multi-Scale Analysis, HSE University, 101000 Moscow, Russia; gsmirnov@hse.ru (G.S.S.)</footnote>

<footnote n="2">² Laboratoire SUBATECH (UMR 6457—IMT Atlantique, Nantes Université, CNRS/IN2P3), 44300 Nantes, France</footnote>

<footnote>* Correspondence: kalinich@subatech.in2p3.fr; Tel.: +33-2-51-85-84-80</footnote>
of a homo-octahedral, 1M trioctahedral mica, and has a monoclinic unit cell with the space group C2/m. However, Pawley et al. [9] used the trigonal 3T polytype structure of TAP for studying its volumetric behavior at high pressures and temperatures.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** TAP supercell (4 × 2 × 2) with water content x = 1. (a) View along the *a*-axis; (b) View along the *b*-axis (the crystallographic directions are shown by the arrows). Green octahedra—MgO, blue tetrahedra—SiO, red spheres—O, white spheres—H. VESTA software [12] was used to visualize the atomistic model.</figcaption></figure>

²⁹Si NMR spectroscopic measurements indicate the presence of Si vacancies in the tetrahedral sheet of TAP [13,14]. Each such vacancy results in the formation of one additional Mg–O–H and three additional Si–O–H groups in the tetrahedral sheet (Figure 2). Thus, each Q²-type Si site (see Figure S1 of the Supplementary Materials) contains silanol groups that donate hydrogen bonds (H-bonds) to the interlayer H₂O molecules. These results suggest that the formation of TAP can involve a defect mechanism that allows favorable H-bonding interaction between the interlayer H₂O molecules and the normally hydrophobic siloxane oxygens of the talc-like layer [14,15]. The proportion of the Q²-type Si sites inferred from the NMR data is estimated to be around 10% [14], corresponding to one silanol for each pair of six-membered siloxane rings of the talc tetrahedral sheet. (Here we are using the common notation for silicon-oxygen tetrahedra, Q'', where the superscript shows the number of other silicon-oxygen tetrahedra attached to the silicon tetrahedron under study [16,17]).

Given the uncertainties in the structure and composition of TAP, computational atom-
istic modeling can be used as a powerful tool to clarify this picture. Classical MD simula-
tions using semi-empirical force fields to describe interatomic interactions within model
systems have made a significant contribution over the last 10–15 years to the detailed
understanding of the structure and properties of clays and clay-related minerals, other
complex nanostructured and nanoporous materials, and their interaction with water and
aqueous solutions [18–23]. ClayFF [18,22] has emerged as one of the most popular and
widely used force fields and has already been thoroughly tested in the atomistic simulations
of such systems [22–24]. Its recent modification, ClayFF-MOH, explicitly takes into account
Metal–O–H (M–O–H) angular bending motions in the mineral structure [25,26], leading
to a better description of the hydroxyl behavior on the edges of mineral particles and
at irregular surfaces [27–29]. At the same time, ab initio methods of atomistic modeling
are now also widely used to study layered minerals and other similar materials [30–40].
Based on a more rigorous quantum chemical foundation of the density functional theory
(DFT) than force-field-based methods, they require, however, orders of magnitude more
computational power for their effective use (e.g., [38]).
<figure><img src="image_6.png" /><figcaption>**Figure 2.** The structural defect in the tetrahedral sheet of TAP. Green—MgO octahedra, blue—SiO₂ tetrahedra, red spheres—O, white spheres—H. (See Section 3.2 for explanation of the labels for each specific type of atoms in the structure.) VESTA software [12] was used to visualize the atomistic model.</figcaption></figure>

Both classical and ab initio methods have been already successfully applied to study the structure and properties of talc, its interaction with water, and an idealized model of TAP, which did not include any structural defects in the talc-like layers [30,39,41,42]. Here, we used classical molecular dynamics (CMD) and ab initio molecular dynamics (AIMD) simulations to study the structure and properties of a more realistic model of TAP that includes structural defects. All CMD simulations were performed both for original and modified versions of the ClayFF force field. Comparing these CMD results with the results of the AIMD simulations and available experimental data allowed us to make a judgment about the reliability and ranges of applicability for each of the computational approaches. All simulations reported here were performed at ambient conditions with the objective of testing the new models and new version of the force field. CMD and AIMD simulations at high temperatures and pressures using these models and analysis of the potential geochemical and geophysical implications of these results will be reported in a separate publication.

**2. Models and Methods**

2.1. *Construction of TAP Models*

The 1M polytype of the phlogopite-type stacking TAP with water content x = 1 was used in this paper as it is estimated to be most stable under ambient conditions [1,41]. The ideal model of TAP was based on experimental single-crystal X-ray diffraction data, which provided unit cell parameters of *a* = 5.323(1) Å, *b* = 9.203(1) Å, *c* = 10.216(1) Å, and *β* = 99.98° with space group C2/m [11]. The supercell containing 6144 atoms (8 × 4 × 4 unit cells along *a*, *b*, and *c* vectors, respectively) was used in the CMD simulations. For the AIMD calculations, smaller supercells with only 768 atoms (4 × 2 × 2 unit cells) were used.

The smaller supercell with 768 atoms was used as the basis for constructing the TAP model with structural defects. First, a random Si atom in a tetrahedral sheet was deleted. Then, four H atoms were added at the site of the Si vacancy to protonate the oxygen atoms that used to coordinate with that Si, forming four hydroxyl groups in its place (Figure 2). To assure a random and uniform distribution of defects, further Si vacancy locations were selected so that the distances between defects in the structure were about 12° Å. The procedure was repeated for all tetrahedral sheets in the supercell, resulting in the creation of one Si defect for every 32 Si atoms in the crystal structure. Such arrangement of defects
closely reproduced the experimentally determined concentration of defects with ~10% of Q²-type Si sites [13,14].

The resulting supercell with structural defects contained 780 atoms and was used in all AIMD calculations. This small supercell was duplicated in all three dimensions to produce a larger supercell with 6240 atoms for CMD calculations, which was equivalent to 8 × 4 × 4 unit cells.

## 2.2. Classical MD Simulations

The original version of the ClayFF force field, ClayFF-orig [18], and its more recent modification, ClayFF-MOH [22], were used to describe the interatomic interactions in two series of simulations in order to make a detailed performance comparison between both versions, and to compare them both with the results of the AIMD simulations described below. The values of the additional parameters for the M–O–H angle bending terms of the ClayFF-MOH version were recently determined by fitting the structural and spectroscopic results of CMD simulations to the results of AIMD calculations for brucite (Mg(OH)₂), gibbsite (Al(OH)₃), and kaolinite (Al₂Si₂O₅(OH)₄) using a simple harmonic functional form [22,25,26]:

$$E_{\text{MOH}} = k(\theta - \theta_0)^2 \qquad (1)$$

where $θ_0$ represents the equilibrium bond angle for a three-body interaction, and k is the stiffness coefficient. The following parameters for Mg–O–H and Si–O–H angles were used in our calculations [22,25,26]: $θ_{0,MgOH} = 110°$, $k_{MgOH} = 6 \text{ kcal} \cdot \text{mol}^{-1} \cdot \text{rad}^{-2}$, $θ_{0,SiOH} = 100°$, and $k_{SiOH} = 15 \text{ kcal} \cdot \text{mol}^{-1} \cdot \text{rad}^{-2}$. In addition, the original harmonic hydroxyl bond stretching terms were replaced here with a more accurate Morse potential [22,42]. Water molecules were described by the flexible SPC/E model [43].

The LAMMPS (7 January 2022 version, Sandia National Laboratories, Albuquerque, NM, USA) software package [44,45] was used to perform all CMD simulations. The cutoff radius for calculating the short-range Lennard–Jones interatomic interactions was 12.5 Å, and the Particle–Particle–Particle–Mesh method was used for the long-range electrostatic interactions [46]. The standard Lorentz–Berthelot combining rules were used to calculate the Lennard–Jones parameters for different atom types [22]. The classical Newtonian equations of atomic motion were numerically integrated with a timestep of 1.0 fs using the velocity Verlet algorithm [46]. The developed atomistic models of TAP were initially equilibrated for 1 ns at 300 K and 1 bar using the Nosé–Hoover thermobarostat [46]. During this equilibration in the NPT statistical ensemble (constant number of particles, pressure, and temperature), no symmetry constraints were imposed on the crystal structures, and all cell parameters were allowed to vary. After the NPT equilibration, the equilibrium simulation run was performed for another 1 ns in the NVT statistical ensemble (constant number of particles, volume, and temperature) using the Nosé–Hoover thermostat [46]. The collected equilibrium dynamic trajectories of the atoms were then used for further statistical analysis.

## 2.3. Ab Initio MD Simulations

The DFT and AIMD calculations were performed using the Gaussian and plane wave basis approach, as implemented in the CP2K simulation software package (version 2022.1, T.D.Kühne et al., EU) [47]. Göedecker–Teter–Hütter (GTH) pseudopotentials [48–50] were used for Mg, Si, O, and H atoms including 10, 4, 6, and 1 valence electrons, respectively. Double-zeta valence polarized (DZVP MOLOPT) basis sets [51] were used for all calculations along with an auxiliary plane wave with the 600 Ry energy cutoff. The generalized gradient approximation (GGA) parametrized by Perdew et al. [52] was used for the exchange-correlation terms with Grimme D3 dispersion correction [53] without the C₉ term. Both the idealized (talc-based) and defect-containing small supercells of TAP were optimized using the Broyden–Fletcher–Goldfarb–Shanno method [54] before starting the AIMD simulations. The nuclear dynamics were treated within the Born–Oppenheimer approximation and the convergence criterion on forces was chosen to be 10⁻⁶ a.u. All AIMD
calculations were performed in the NVT ensemble at 300 K with a 0.5 fs timestep using the Nosé-Hoover scheme [47]. Each system was pre-equilibrated for 2 ps before the 10 ps production run was performed for further statistical analysis of the resulting properties.

Periodic boundary conditions [46,47] were applied in all CMD and AIMD calculations and no symmetry constraints were imposed on the simulated structures.

### 3. Results and Discussion

#### 3.1. Crystallographic Parameters

The equilibrium crystallographic unit cell parameters of the idealized and defective TAP structures were obtained using CMD calculations at 300 K and from the zero temperature cell optimization with DFT (Table 1), as described in the previous section. Six models were considered in total, all of which were in qualitative agreement with the experimental data [7,10,11,13,55]. Previously, only the ideal talc-based TAP structure was investigated using CMD simulations with the ClayFF-orig force field [41] and DFT calculations without dispersion corrections [30].

<table><caption>**Table 1.** Crystallographic unit cell parameters of the TAP models: Simulated results and experimental data.</caption><thead><tr><th>Source</th><th>a, Å</th><th>b, Å</th><th>c, Å</th><th>V, Å<sup>3</sup></th><th>α, °</th><th>β, °</th><th>γ, °</th></tr></thead><tbody><tr><td>ClayFF-orig (Ideal)</td><td>5.288</td><td>9.169</td><td>9.966</td><td>475.2</td><td>90.03</td><td>100.38</td><td>89.99</td></tr><tr><td>ClayFF-MOH (Ideal)</td><td>5.291</td><td>9.175</td><td>9.969</td><td>476.0</td><td>90.01</td><td>100.36</td><td>89.98</td></tr><tr><td>DFT (PBE + D3) (Ideal)</td><td>5.357</td><td>9.251</td><td>9.961</td><td>486.2</td><td>90.00</td><td>99.98</td><td>90.00</td></tr><tr><td>ClayFF-orig (Defect)</td><td>5.308</td><td>9.206</td><td>9.945</td><td>478.0</td><td>89.99</td><td>100.38</td><td>90.00</td></tr><tr><td>ClayFF-MOH (Defect)</td><td>5.316</td><td>9.221</td><td>9.939</td><td>479.3</td><td>89.97</td><td>100.31</td><td>90.00</td></tr><tr><td>DFT (PBE + D3) (Defect)</td><td>5.376</td><td>9.276</td><td>9.943</td><td>488.3</td><td>90.00</td><td>99.98</td><td>90.00</td></tr><tr><td>ClayFF (Ideal) [41]</td><td>5.288</td><td>9.17</td><td>9.99</td><td>477.1</td><td>90</td><td>100</td><td>90</td></tr><tr><td>DFT [30]</td><td>5.273</td><td>9.137</td><td>9.397</td><td>453.2</td><td>90</td><td>99.99</td><td>90</td></tr><tr><td>XRD [10]</td><td>5.316</td><td>9.191</td><td>10.118</td><td>486.5</td><td>90</td><td>100.25</td><td>90</td></tr><tr><td>XRD [7]</td><td>5.293</td><td>9.194</td><td>10.044</td><td>486.0</td><td>90</td><td>96.11</td><td>90</td></tr><tr><td>XRD [11]</td><td>5.323</td><td>9.203</td><td>10.216</td><td>492.9</td><td>90</td><td>99.98</td><td>90</td></tr><tr><td>XRD [13]</td><td>5.3297</td><td>9.205</td><td>10.202</td><td>492.8</td><td>90</td><td>100.82</td><td>90</td></tr><tr><td>XRD [55]</td><td>5.3007</td><td>9.186</td><td>10.185</td><td>488.5</td><td>90</td><td>99.97</td><td>90</td></tr></tbody></table>

Among the classical models, the best agreement was achieved with the defective structure and the ClayFF-MOH force field; the error of the unit cell volume did not exceed 3% (Table 1). The results for the idealized structure with ClayFF-orig virtually coincided with the previous calculations [41].

The zero temperature DFT calculations gave the most accurate value of the unit cell volume compared to the experimental data, but the cell vectors and shape of the unit cell were not quite correct. The cell was squeezed along the c-axis and stretched along the a- and b-axes. The difference between the ideal and defect models was rather small, probably due to the symmetry constraints during the DFT optimization procedure.

The lattice parameters of the earlier DFT calculations [30] significantly deviated from the experimental data. However, those results were obtained without a dispersion correction, which is especially important in systems containing hydroxyls and/or H₂O molecules in the structure [25,56].
## 3.2. Atom Positions and Interatomic Distances

All calculated average distances between metal atoms (Mg, Si) and various oxygen atoms in the structure (O<sub>a</sub>—apical oxygen atom; O<sub>b</sub>—basal oxygen; O<sub>h</sub>—oxygen atom of hydroxyl; O<sub>w</sub>—oxygen atom of H<sub>2</sub>O) were in fairly good agreement with the experimental results [11,55] (Table 2). The maximum deviation of the calculated values from the experimental data was below 6%. However, the AIMD calculations usually better reproduced the experimental bond distances compared to the CMD simulations.

<table><caption>**Table 2.** Average distances between some atoms and the interlayer thickness of TAP: Simulated results and experimental data.</caption><thead><tr><th>Source</th><th>Si–O<sub>b</sub>, Å</th><th>Si–O<sub>a</sub>, Å</th><th>Mg–O<sub>a</sub>, Å</th><th>Mg–O<sub>h</sub>, Å</th><th>O<sub>w</sub>–O<sub>b</sub>, Å</th><th>Interlayer Thickness, Å</th></tr></thead><tbody><tr><td>ClayFF-orig (Ideal)</td><td>1.576</td><td>1.527</td><td>2.075</td><td>2.162</td><td>3.158</td><td>3.498</td></tr><tr><td>ClayFF-MOH (Ideal)</td><td>1.577</td><td>1.526</td><td>2.076</td><td>2.153</td><td>3.161</td><td>3.507</td></tr><tr><td>DFT (PBE + D3) (Ideal)</td><td>1.622</td><td>1.609</td><td>2.076</td><td>2.050</td><td>3.151</td><td>3.339</td></tr><tr><td>ClayFF-orig (Defect)</td><td>1.583 [Si<sub>s</sub>–O<sub>bs</sub>] * 1.576</td><td>1.520 [Si<sub>s</sub>] 1.526</td><td>2.081</td><td>2.093 [O<sub>hs</sub>] 2.163</td><td>2.935 [O<sub>bs</sub>] 3.171</td><td>3.502</td></tr><tr><td>ClayFF-MOH (Defect)</td><td>1.579 [Si<sub>s</sub>–O<sub>bs</sub>] 1.577</td><td>1.516 [Si<sub>s</sub>] 1.525</td><td>2.082</td><td>2.095 [O<sub>hs</sub>] 2.156</td><td>2.708 [O<sub>bs</sub>] 3.199</td><td>3.519</td></tr><tr><td>DFT (PBE + D3) (Defect)</td><td>1.631 [Si<sub>s</sub>–O<sub>bs</sub>] 1.623</td><td>1.613 [Si<sub>s</sub>] 1.607</td><td>2.086</td><td>1.999 [O<sub>hs</sub>] 2.053</td><td>2.778 [O<sub>bs</sub>] 3.179</td><td>3.348</td></tr><tr><td>XRD [11]</td><td>1.633</td><td>1.607</td><td>2.089</td><td>2.064</td><td>3.192</td><td>3.526</td></tr><tr><td>XRD [55]</td><td>1.626</td><td>1.617</td><td>2.076</td><td>2.056</td><td>3.190</td><td>3.525</td></tr></tbody></table>

Atomistic modeling of the TAP structure with defects made it possible to obtain the distances between atoms, at least one of which belonged to a structural defect. The average distance between Si<sub>s</sub> and O<sub>bs</sub> (Si and basal O belonging to silanol groups) was slightly larger than the Si–O<sub>b</sub> distance in all simulations. The same result was also obtained for the Si–O<sub>a</sub> pair in the AIMD calculations (Table 2), although the CMD simulations gave a slightly smaller Si<sub>s</sub>–O<sub>a</sub> length.

Slightly smaller interatomic distances were also observed in all simulations for the Mg–O<sub>hd</sub> (hydroxyl oxygen inside the defect) and O<sub>w</sub>–O<sub>bs</sub> pairs compared to the Mg–O<sub>h</sub> and O<sub>w</sub>–O<sub>b</sub> couples, respectively. Such differences could be explained by the formation of H-bonds donated by H<sub>hs</sub> (H of silanol groups) to O<sub>w</sub> and by H<sub>hd</sub> (hydroxyl hydrogen inside the defect) to O<sub>bs</sub>.

The experimental interlayer thickness of TAP was poorly reproduced with the AIMD simulations; it was about 0.2 Å lower in both cases. A similar discrepancy was also observed in the DFT simulations of talc without dispersion corrections [31]. The best agreement with the experimental data for TAP was observed using the defect model and the ClayFF-MOH force field.

## 3.3. Radial Distribution Functions

Radial distribution functions (RDFs) are especially important for understanding the structural arrangement and ordering of the H<sub>2</sub>O molecules in the TAP interlayers. They were calculated for different H–O pairs existing in our models. The first H<sub>w</sub>–O<sub>w</sub> peaks from the AIMD calculations were located at 4.7 Å. The results were almost the same for the ideal and defect models (Figure 3a; see also Figure 2), probably due to interlayer shrinkage leading to less mobile H<sub>2</sub>O molecules compared to the experimental data (Table 2).

<footnote>* The atoms belonging to the structural defect are additionally indicated in square brackets.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Radial distribution functions. (a) H<sub>w</sub>–O<sub>w</sub> pair; (b) H<sub>w</sub>–O<sub>b</sub> pair. See Figure 2 for notations.</figcaption></figure>

The first peak for the ideal model was slightly shifted to the left in the CMD simula-
tions, but the AIMD-like position was obtained using the defect model with the ClayFF-
MOH force field. This was due to the formation of stronger H-bonds between the hydrogen
atom of the silanol groups (H<sub>hs</sub>) and O<sub>w</sub> of the water molecules (see also Figure 2). The
hydroxyl groups in the structural defects enhanced the hydrophilicity, exhibiting stronger
attraction towards the nearest H<sub>2</sub>O molecules approaching the tetrahedral layer. Conse-
quently, the nearest water molecules were slightly displaced away from other interlayer
water molecules. This effect was not observed with the ClayFF-orig force field.

The position of the first H<sub>w</sub>–O<sub>b</sub> RDF peak for the ClayFF-MOH model with defects was also in good agreement with the AIMD calculations (Figure 3b; see also Figure 2). However, all RDF peaks from the CMD simulations were broader and the average H<sub>w</sub>–O<sub>b</sub> distance from the AIMD calculations was smaller than the classical values (3.3 Å vs. 3.4–3.6 Å, respectively). This was also consistent with the smaller interlayer space thickness resulting from the AIMD calculations (see Table 2).

The H<sub>h</sub>–O<sub>w</sub> RDFs (Figure 4a; see also Figure 2) showed different widths for the first peak. The first peak for the ClayFF-MOH model with defects was wider compared to other classical models and closer to the AIMD results. This could be explained by the greater mobility of the H<sub>2</sub>O molecules in the six-membered tetrahedral rings [11], which were free from short-range interaction due to the structural defects.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Radial distribution functions. (a) H<sub>h</sub>–O<sub>w</sub> pair; (b) H<sub>hs</sub>–O<sub>w</sub> pair (defect model). See Figure 2 for notations.</figcaption></figure>

This suggestion was also confirmed by the H<sub>hs</sub>-O<sub>w</sub> RDFs for the defective structures
(Figure 4b; see also Figure 2), where the first peak locations resulting from the ClayFF-MOH
and AIMD simulations were almost the same (1.8 and 1.7 Å, respectively). Such short
donor-acceptor distances, especially in the case of the AIMD simulations, indicated the
formation of particularly strong H<sub>hs</sub>...O<sub>w</sub> H-bonds compared to the CMD calculations with
ClayFF-orig.

Fairly good agreement with the ClayFF-MOH and AIMD simulation results was also observed for the H<sub>hd</sub>-O<sub>bs</sub> RDFs (Figure 5; see also Figure 2). The average distance between H<sub>hd</sub> and O<sub>bs</sub> from the simulations with the ClayFF-orig model was smaller compared to the results of the ClayFF-MOH and AIMD simulations. The peak intensity for this RDF from the ClayFF-MOH and AIMD simulations was also higher than for the ClayFF-orig results. This possibly indicated competition between three O<sub>bs</sub> atoms for H-bond formation with one of the H<sub>hd</sub> atoms inside the structural defect for the ClayFF-MOH and AIMD models, leading to the formation of so-called bifurcated H-bonding [57]. For the ClayFF-orig model, the formation of a single stronger H-bond between one of H<sub>hd</sub> and O<sub>bs</sub> was most likely.

<figure><img src="image_9.png" /><figcaption>Figure 5. Radial distribution functions for the H<sub>hd</sub>-O<sub>bs</sub> pair (defect model).</figcaption></figure>

3.4. *Interlayer Atomic Density Distributions*

To better understand the TAP interlayer structure, 2-dimensional contour maps of
the time-averaged atomic density distributions were calculated from the CMD and AIMD
simulations (Figures 6 and 7).

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Time-averaged atomic density distributions of the interlayer space of TAP models with structural defects from CMD. (a) ClayFF-orig; (b) ClayFF-MOH. Green contours—O<sub>b</sub>; black contours—H<sub>w</sub>; red contours—O<sub>w</sub>; brown contours—H<sub>hs</sub>; magenta contours—O<sub>bs</sub>; blue contours—H<sub>hd</sub>.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Time-averaged atomic density distributions of the interlayer space of TAP from AIMD. (a) Ideal model; (b) Defect model. Green contours—O<sub>b</sub>; black contours—H<sub>w</sub>; red contours—O<sub>w</sub>; brown contours—H<sub>hs</sub>; magenta contours—O<sub>bs</sub>; blue contours—H<sub>hd</sub>.</figcaption></figure>

For the ideal TAP model, there were no differences between the two versions of ClayFF and AIMD results. In all cases, the H₂O molecules were located above the structural hydroxyls of each six-membered siloxane ring [41] of the talc tetrahedral layers (see Figure 7a). Also, the H₂O molecules were coordinated by the O<sub>b</sub> atoms of the tetrahedral layer [41]. According to previous DFT results [30], the H<sub>w</sub> atom can form multi-furcated H-bonds with O<sub>b</sub> atoms.

Water molecules located far from the structural defects had similar behaviors to the H₂O molecules in the ideal model for both ClayFF versions (Figure 6). However, the H₂O molecules near the defects behaved differently due to their stronger interaction with the silanol groups of the defects. The CMD results with the ClayFF-MOH model demonstrated an ordered pattern of H₂O molecule arrangement near the defects, which was due to their acceptance of strong H-bonds from the silanol groups of the defects (Figure 7b).

### 3.5. *Hydrogen Bonding Structure and Dynamics in the TAP Interlayers*

A common geometrical definition was used to determine whether a H-bond (HB) existed between a donor and acceptor [58,59]: R<sub>OdOa</sub> ≤ 3.5 Å, R<sub>OaHd</sub> ≤ 2.45 Å, and φ<sub>HdOdOa</sub> ≤ 30°, where *d* is a donor and *a* is an acceptor. The average lifetime of H-bonds (τ<sub>HB</sub>) was estimated by integrating the so-called continuous time functions of H-bonds [58,59]. The average number of H-bonds (n<sub>HB</sub>) per donor was also calculated (Table 4; see also Figures S2 and S3 of the Supplementary Materials).

<table><caption>**Table 3.** The average number of H-bonds per donor atom, n<sub>HB</sub>, and the average lifetime of H-bonds, τ<sub>HB</sub>, for different TAP models.</caption><thead><tr><th rowspan="2">Donor-Acceptor Pair</th><th colspan="3">n<sub>HB</sub></th><th colspan="3">τ<sub>HB</sub>, fs</th></tr><tr><th>ClayFF-orig</th><th>ClayFF-MOH</th><th>AIMD</th><th>ClayFF-orig</th><th>ClayFF-MOH</th><th>AIMD</th></tr></thead><tbody><tr><td>H<sub>w</sub>···O<sub>b</sub></td><td>0.95/0.88 *</td><td>0.90/0.70</td><td>0.75/0.73</td><td>21/23</td><td>20/22</td><td>19/21</td></tr><tr><td>H<sub>h</sub>···O<sub>w</sub></td><td>0.19/0.19</td><td>0.19/0.25</td><td>0.10/0.13</td><td>76/88</td><td>77/143</td><td>61/96</td></tr><tr><td>H<sub>h</sub>···O<sub>b</sub></td><td>0.02/0.01</td><td>-</td><td>-</td><td>6/6</td><td>-</td><td>-</td></tr><tr><td>H<sub>hs</sub>···O<sub>w</sub></td><td>0.30</td><td>0.86</td><td>0.77</td><td>29</td><td>249</td><td>1453</td></tr></tbody></table>
<table><caption>**Table 4.** *Cont.*</caption><thead><tr><th rowspan="2">Donor-Acceptor Pair</th><th colspan="3">n<sub>HB</sub></th><th colspan="3">τ<sub>HB</sub>, fs</th></tr><tr><th>ClayFF-orig</th><th>ClayFF-MOH</th><th>AIMD</th><th>ClayFF-orig</th><th>ClayFF-MOH</th><th>AIMD</th></tr></thead><tbody><tr><td>H<sub>w</sub>···O<sub>bs</sub></td><td>0.14</td><td>0.29</td><td>-</td><td>74</td><td>210</td><td>-</td></tr><tr><td>H<sub>hs</sub>···O<sub>bs</sub></td><td>0.20</td><td>-</td><td>-</td><td>18</td><td>-</td><td>-</td></tr><tr><td>H<sub>hd</sub>···O<sub>bs</sub></td><td>0.48</td><td>0.25</td><td>0.49</td><td>14</td><td>8</td><td>8</td></tr></tbody></table>

### 3.5.1. Ideal Model

Three possible donor-acceptor pairs in the ideal TAP model were possible: H<sub>w</sub>···O<sub>b</sub> (H<sub>2</sub>O—basal oxygen atoms), H<sub>h</sub>···O<sub>w</sub> (structural hydroxyl—H<sub>2</sub>O), and H<sub>h</sub>···O<sub>b</sub> (structural hydroxyl—basal oxygen atoms). A previous theoretical study [30] revealed the formation of only two donor-acceptor pairs, H<sub>w</sub>···O<sub>b</sub> and H<sub>h</sub>···O<sub>w</sub>. Our CMD results with the ClayFF-MOH force field and AIMD results confirmed the presence of only these two donor-acceptor pairs (Table 4). Weak H<sub>h</sub>···O<sub>b</sub> H-bonds occurred only in the CMD simulation with ClayFF-orig, which was due to the absence of the Mg-O-H angle bending term in this earlier model. The lifetimes of these H-bonds were generally short and did not exceed 80 fs. However, the H<sub>h</sub>···O<sub>w</sub> lifetime was the longest compared to the lifetime of H<sub>w</sub>···O<sub>b</sub> H-bonds. The average number of H-bonds was similar in all CMD simulations, while the AIMD results predicted lower values.

The short lifetimes of H-bonds indicated that the H<sub>2</sub>O molecules had significant rotational mobility. This suggestion was confirmed by the second-order orientational time correlation function [60,61], which was calculated using the orientation of the unit vector along the O<sub>w</sub>-H<sub>w</sub> bond in the H<sub>2</sub>O molecules (Figure 8). Thus, the interaction between structural hydroxyls and H<sub>2</sub>O molecules played a very important role in the behavior of the interlayer H<sub>2</sub>O molecules in the ideal model.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Autocorrelation function of H<sub>2</sub>O molecule reorientation.</figcaption></figure>

### 3.5.2. Defect Model

The observed behavior of the H<sub>2</sub>O molecules far from the TAP structural defects was the same as in the ideal model. The H<sub>w</sub>···O<sub>b</sub> H-bonding lifetime did not change significantly (Table 4). However, the number of H<sub>w</sub>···O<sub>b</sub> H-bonds was smaller in all simulations. On the contrary, an increase in the average number of H<sub>h</sub>···O<sub>w</sub> H-bonds was observed for the ClayFF-MOH and AIMD simulations. The lifetime of H<sub>h</sub>···O<sub>w</sub> H-bonds was also increased, and the longest lifetime was observed for the ClayFF-MOH force field (see Figure S2b of the Supplementary Material). This strengthening of H-bonds was also reflected in the orientational relaxation of the H<sub>2</sub>O molecules. The orientational relaxation time was the

<footnote>* The first value is for the ideal model; the second value is for the model with defects.</footnote>
largest for the defect model with the ClayFF-MOH force field (Figure 8). Similar to the ideal case, $H_h...O_b$ H-bonds were formed only with the ClayFF-orig force field.

In the new TAP model with defects, several new donor-acceptor H-bonding pairs were possible in addition to the main $H_w...O_b$ and $H_h...O_w$ H-bonds. The most important was the $H_{hs}...O_w$ pair. It had the longest lifetime, especially as observed in the AIMD simulations (Table 4 and Figure S3 of the Supplementary Materials). However, the simulation results with the ClayFF-orig force field led to only very weak H-bonds for this donor-acceptor pair. This was due to the possibility of frequent and completely unrestricted reorientation of the silanol hydroxyl groups. The average number of $H_{hs}...O_w$ H-bonds was also the smallest for the ClayFF-orig case, while fairly close results were observed for the ClayFF-MOH and AIMD simulations.

Another important H-bonding pair was $H_w...O_{bs}$, which was due to H-bond donation by the $H_2O$ molecules to the oxygen atoms of the silanols. The CMD simulations demonstrated the existence of such H-bonds, but they were not observed during the AIMD runs. The stronger $H_{hs}...O_w$ H-bonds probably prevented the formation of such stable pairs in the AIMD simulations. The formation of weak H-bonds donated by the structural hydroxyls to the oxygen atoms of the silanol groups ($H_{hd}...O_{bs}$ pairs) was also possible inside the structural defects. Thus, the interaction of the silanol groups with the interlayer $H_2O$ molecules played even a greater role in the defect TAP model.

### 3.6. Vibrational Properties

The power spectra of atomic motion or vibrational density of state (VDOS) were obtained by Fourier transform of the velocity autocorrelation function of the respective atoms (see, e.g., [23]).

The total VDOS (Figure S4 of the Supplementary Materials) was decomposed into several atomic and molecular contributions: $H_2O$ molecules (Figure S5 of the Supplementary Materials), $H_h$ atoms (Figure S6 of the Supplementary Materials), and $H_{hs}$ atoms (Figure 9). According to Raman spectroscopic data in the high-frequency region, the peak of the $O_h-H_h$ stretching mode was located at 3622 $cm^{-1}$, and there were two peaks reflecting the symmetric and asymmetric stretching of the $O_w-H_w$ bonds within the water molecules at 3593 $cm^{-1}$ and 3668 $cm^{-1}$ [1].

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Power spectra of hydrogen atoms ($H_{hs}$) of the silanol groups of the defect TAP model.</figcaption></figure>

The stretching modes of the $H_2O$ molecules obtained from the CMD simulations were located at slightly higher frequencies due to the fact that only a simple harmonic function for the intramolecular vibrations was used in the SPC/E water model (Table 5). A more accurate and elaborate $H_2O$ intramolecular potential would be necessary to better reproduce the positions and widths of these peaks (e.g., [62]). The $O_w-H_w$ stretching modes
from the AIMD calculations showed better agreement with the experimental data; however, the resolution was not really sufficient to clearly distinguish the two peaks.

<table><caption>**Table 5.** Calculated and experimental vibrational peak positions in the high-frequency region for TAP models at ambient conditions.</caption><thead><tr><th rowspan="3">Modes</th><th colspan="7">Wavenumber, cm<sup>-1</sup></th></tr><tr><th colspan="3">Ideal Model</th><th colspan="3">Defect Model</th><th rowspan="2">Experiment [1]</th></tr><tr><th>ClayFF-orig</th><th>ClayFF-MOH</th><th>AIMD</th><th>ClayFF-orig</th><th>ClayFF-MOH</th><th>AIMD</th></tr></thead><tbody><tr><td>O–H<sub>sym</sub> stretching of interlayer H<sub>2</sub>O</td><td>3765</td><td>3765</td><td>3727</td><td>3761</td><td>3748</td><td>3663</td><td>3593</td></tr><tr><td>O–H<sub>asym</sub> stretching of interlayer H<sub>2</sub>O</td><td>3823</td><td>3823</td><td>3771</td><td>3827</td><td>3820</td><td>3728</td><td>3668</td></tr><tr><td>O–H stretching of structural hydroxyls</td><td>3569</td><td>3569</td><td>3711</td><td>3573</td><td>3586</td><td>3703</td><td>3622</td></tr><tr><td>O–H stretching of silanols</td><td>-</td><td>-</td><td>-</td><td>3577</td><td>3431</td><td>3194</td><td>-</td></tr></tbody></table>

The AIMD peak positions reflecting the stretching mode of structural hydroxyls (H<sub>h</sub> atoms) were blue-shifted compared to the experimental data and the results of the CMD simulations (Table 5). The experimentally observed peak at 3622 cm<sup>-1</sup> was better reproduced for the TAP model with defects using the ClayFF-MOH classical force field (see Figure S6 of the Supplementary Materials).

In the high-frequency region, there were also peaks reflecting the O–H stretching mode of the hydroxyls of the silanol groups (H<sub>hs</sub> atoms) (Figure 9). The peak position obtained from the AIMD simulations was located at 3194 cm<sup>-1</sup>, while the CMD results showed peaks located at 3577 and 3431 cm<sup>-1</sup> for ClayFF-orig and ClayFF-MOH, respectively (Table 5). This was due to stronger H-bond formation between the hydroxyls of the silanol groups and the interlayer H<sub>2</sub>O molecules (see Table 4).

In all calculations, the vibrational peaks of the silanol groups were located at lower frequencies than similar peaks at the edges of the talc crystals [63]. They were also red-shifted compared to the high-pressure experimental data of Pawley and Welch [64], in which these vibrational modes were detected at 3587 cm<sup>-1</sup>. However, these modes are sensitive to pressure and the strength of interlayer H-bonds, thus additional studies are required for a more detailed quantitative comparison with the experimental data.

The Raman spectrum of TAP in the lower frequency region [1] shows several vibrational modes assigned to Mg–O–H and Si–O–Si bending, OH transition, OH libration, Si–O–Si symmetric stretching, and Si–O stretching. The VDOS calculated from our CMD and AIMD simulations showed good agreement with some experimental peaks in this lower frequency region. Generally, we were able to identify the OH librations in the calculated spectra.

The calculated partial VDOS for H<sub>h</sub> from the CMD simulations with ClayFF-orig in models with and without defects showed that the associated peak was located at a lower frequency compared to the ClayFF-MOH and AIMD models (Figure S6 of the Supplementary Materials). This was a clear indication that the addition of the constraining Mg–O–H angular bending term in the ClayFF-MOH model induced better agreement with the AIMD results. Better agreement with the ClayFF-MOH and AIMD results was also observed for the Si–O–H bending vibrations of the silanol groups of the defects in the lower frequency region of spectra, reflecting librational motions (Figure 9). Such librational modes of H<sub>2</sub>O were located at ca. 310–320 cm<sup>-1</sup> in the CMD simulations and at ca. 375–385 cm<sup>-1</sup> in the AIMD simulations (Figure S5 of the Supplementary Materials). The results were consistent with the Raman spectroscopic data [1,65].
In general, due to the intrinsic uncertainty of these kinds of classical model, a comparison of absolute values of the calculated vibrational frequencies with the experimental data are less meaningful than a comparison of trends in these vibrational modes with temperature, pressure, and composition (e.g., [22,23,29]). These trends deserve a much more detailed analysis and discussion and will be the subject of a separate publication.

### 3.7. Elastic Properties

CMD simulations with the ClayFF force field have been successfully used to calculate the elastic and thermal properties of clays and clay-related materials [22,66,67], and it has been recently demonstrated that the modified ClayFF-MOH version allows for better reproduction of elastic constants for a number of minerals [28,29].

A series of special equilibrium CMD runs was performed here to calculate the elastic constants of different TAP models. First, the stress tensor components were calculated as the sum of the kinetic and virial terms for the equilibrium supercell. Then, negative and positive supercell deformations of 1.0% were applied to the supercell in six independent directions, and the stress tensor components were calculated for the specified deformations. Such amplitudes of deformation are suitable for talc-like materials [37]. The CMD calculations of the initial and deformed TAP supercells were carried out in the NVT statistical ensemble for 0.5 ns. The elastic constants of the crystal were then numerically determined according to the generalized Hooke's law. The bulk modulus ($K_H$) and shear modulus ($G_H$) of TAP were also calculated using the Voigt-Reuss-Hill approximation [68].

It is important to emphasize that all components of the elastic tensor were calculated completely independently without any symmetry constraints imposed on the crystal structure (Table 6). They were compared with previous DFT calculations for monoclinic talc (C2/c) [37,69]. The off-diagonal components were less than 10 GPa in all cases. The TAP and talc elastic constants were very close to each other, except for the $C_{33}$ and $C_{55}$ values. The $C_{33}$ constant is responsible for the elasticity along the c-axis. Smaller values for TAP are, obviously, due to the presence of H$_2$O molecules in the interlayers, which leads to a softer behavior along the c-axis.

<table><caption>**Table 6.** Elastic properties of TAP models and talc (monoclinic) under ambient conditions in Voigt notation.</caption><thead><tr><th>Model (Source)</th><th>C<sub>11</sub></th><th>C<sub>12</sub></th><th>C<sub>13</sub></th><th>C<sub>15</sub></th><th>C<sub>22</sub></th><th>C<sub>23</sub></th><th>C<sub>25</sub></th><th>C<sub>33</sub><br>GPa</th><th>C<sub>35</sub></th><th>C<sub>44</sub></th><th>C<sub>46</sub></th><th>C<sub>55</sub></th><th>C<sub>66</sub></th><th>K<sub>H</sub></th><th>G<sub>H</sub></th></tr></thead><tbody><tr><td>Ideal TAP<br>(ClayFF-orig)</td><td>274.1</td><td>120.0</td><td>3.5</td><td>0.8</td><td>273.9</td><td>15.4</td><td>-5.3</td><td>27.1</td><td>-0.4</td><td>2.5</td><td>0.8</td><td>11.7</td><td>78.9</td><td>39.9</td><td>16.9</td></tr><tr><td>Ideal TAP<br>(ClayFF-MOH)</td><td>282.6</td><td>116.7</td><td>12.9</td><td>2.1</td><td>271.1</td><td>7.6</td><td>2.3</td><td>8.8</td><td>2.4</td><td>6.7</td><td>2.2</td><td>7.2</td><td>86.8</td><td>38.9</td><td>19.3</td></tr><tr><td>Defect TAP<br>(ClayFF-orig)</td><td>237.9</td><td>98.6</td><td>10.5</td><td>6.8</td><td>239.9</td><td>15.2</td><td>3.8</td><td>23.1</td><td>4.8</td><td>4.1</td><td>-0.4</td><td>4.1</td><td>63.8</td><td>46.2</td><td>22.6</td></tr><tr><td>Defect TAP<br>(ClayFF-MOH)</td><td>243.6</td><td>92.9</td><td>7.4</td><td>5.3</td><td>240.9</td><td>9.5</td><td>0.5</td><td>10.6</td><td>3.5</td><td>4.9</td><td>-3.7</td><td>1.5</td><td>76.4</td><td>44.4</td><td>23.8</td></tr><tr><td>Talc [69]</td><td>243.2</td><td>96.4</td><td>23.6</td><td>-29.9</td><td>256.6</td><td>10.3</td><td>-15.2</td><td>53.7</td><td>-19.1</td><td>40.1</td><td>-5.5</td><td>54.9</td><td>74.5</td><td>65.0</td><td>54.6</td></tr><tr><td>Talc [37]</td><td>215.8</td><td>68.5</td><td>12.1</td><td>3.2</td><td>218.1</td><td>2.1</td><td>2.6</td><td>63.4</td><td>1.1</td><td>2.9</td><td>0.9</td><td>30.3</td><td>73.8</td><td>60.4</td><td>30.3</td></tr></tbody></table>

The presence of interlayer H$_2$O molecules and the absence of shifts of the TOT layers with respect to each other in all TAP models explained the smaller values of the $C_{55}$ constant, which is responsible for shear elastic behavior. The values of $K_H$ and $G_H$ for the defect TAP model were closer to the talc DFT results. As demonstrated in the previous sections, a stable H-bonding network is formed inside the interlayer space of TAP, which additionally stiffened the defect structure.

### 4. Conclusions

Six atomistic models of the 1M polytype of the 10-Å phase (TAP), Mg$_3$Si$_4$O$_{10}$(OH)$_2 \cdot xH_2O$, were constructed and quantitatively studied using classical and ab initio molecular dynamics.
ics simulations. A water content of $x = 1$ was used in this study, as it was the most stable model under ambient conditions.

The simulation results clearly demonstrated that the inclusion of Si vacancies in the TAP structure and the presence of silanol groups around these structural defects provided a better description of the experimental volumetric data. The advantage of the defect TAP model was also demonstrated by comparing the elastic constants and bulk and shear moduli between the 10-Å phase and monoclinic talc.

Orientational ordering of the hydroxyls of the silanol groups ($O_{bs}-H_{hs}$) was observed in the AIMD simulations, with atomic positions having distinct structural pattern. However, reorientation events were very rare. CMD simulations using the recently modified ClayFF-MOH force field resulted in a similar behavior, while the earlier version of the force field, ClayFF-orig, led to much more mobile and disordered $O_{bs}-H_{hs}$ groups. These structural observations were also supported by the analysis of the vibrational spectra calculated from the results of both the CMD and AIMD simulations. In all cases, the modified ClayFF-MOH version of the force field provided more accurate descriptions that were closer to the DFT and AIMD results than the original version of ClayFF. However, the stretching modes of the $O_{bs}-H_{hs}$ bonds in the CMD simulations were located at somewhat lower frequencies than suggested by the experimental data.

In the idealized talc-based and defectless TAP structure, all types of H-bonds are quite weak, and this facilitates the orientational relaxation of the interlayer H₂O molecules, which is very fast compared to normal liquid water. However, stronger H-bonds were formed between the silanol groups of the structural defects and the interlayer H₂O molecules in the TAP structural model with defects. Thus, a stable H-bonding network could be observed in the interlayers of this TAP model. This phenomenon could certainly play an important role in the behavior of TAP at high pressures and temperatures and could affect the retention and transport of water by this phase in the Earth's upper mantle in subduction zones. A detailed analysis of the CMD and AIMD simulations at high temperatures and pressures for the new model of the 10-Å phase with defects, as well as discussion of the possible geochemical and geophysical implications of these results, will be the subject of a separate publication improving on the earlier simulations of the idealized talc-based and defectless TAP model [41,42].

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/min13081018/s1: Figure S1. Q³- and Q²-type Si sites in the tetrahedral sheet of the TAP; Figure S2. Continuous H-bonding time autocorrelation functions; Figure S3. Continuous H-bonding time autocorrelation functions for the $H_{hs}...O_{w}$ pairs; Figure S4. Total VDOS for the ideal and defect TAP models; Figure S5. Partial VDOS of the interlayer H₂O molecules for the ideal and defect TAP models; Figure S6. Partial VDOS of the $H_{h}$ atoms for the ideal and defect TAP models.

**Author Contributions:** Conceptualization, A.G.K.; methodology, A.G.K. and G.S.S.; investigation, E.V.T.; data curation, E.V.T., G.S.S. and A.G.K.; writing—original draft preparation, E.V.T.; writing—review and editing, E.V.T., G.S.S. and A.G.K.; visualization, E.V.T.; supervision, G.S.S. and A.G.K.; project administration, G.S.S. and A.G.K.; funding acquisition, G.S.S. and A.G.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the HSE University Basic Research Program (E.V.T. and G.S.S.) and ANDRA, Orano, and EDF (A.G.K.) via the industrial chair "Storage and Disposal of Radioactive Waste" at the IMT Atlantique.

**Data Availability Statement:** All data are available from the authors upon request.

**Acknowledgments:** The article was prepared within the framework of the HSE University Basic Research Program and was supported by computational resources of HPC facilities at HSE University [70]. A.G.K. also acknowledges the financial support of the industrial chair "Storage and Disposal of Radioactive Waste" at the IMT Atlantique, funded by ANDRA, Orano, and EDF.

**Conflicts of Interest:** The authors declare no conflict of interest.
References

1. Fumagalli, P.; Stixrude, L.; Poli, S.; Snyder, D. The 10Å phase: A high-pressure expandable sheet silicate stable during subduction of hydrated lithosphere. *Earth Planet. Sci. Lett.* **2001**, *186*, 125–141. [CrossRef]

2. Ohtani, E.; Litasov, K.; Hosoya, T.; Kubo, T.; Kondo, T. Water transport into the deep mantle and formation of a hydrous transition zone. *Phys. Earth Planet Inter.* **2004**, *143–144*, 255–269. [CrossRef]

3. Khisina, N.R.; Wirth, R. Nanoinclusions of high-pressure hydrous silicate, Mg₃Si₄O₁₀(OH)₂·nH₂O (10Å-phase), in mantle olivine: Mechanisms of formation and transformation. *Geochem. Int.* **2008**, *46*, 319–327. [CrossRef]

4. Chollet, M.; Daniel, I.; Koga, K.T.; Petitgirard, S.; Morard, G. Dehydration kinetics of talc and 10Å phase: Consequences for subduction zone seismicity. *Earth Planet Sci. Lett.* **2009**, *284*, 57–64. [CrossRef]

5. Rashchenko, S.V.; Likhacheva, A.Y.; Goryainov, S.V.; Krylov, A.S.; Litasov, K.D. In situ spectroscopic study of water intercalation into talc: New features of 10Å phase formation. *Amer. Mineral.* **2016**, *101*, 431–436. [CrossRef]

6. Yamamoto, K.; Akimoto, S. The system MgO-SiO₂-H₂O at high pressures and temperatures; stability field for hydroxyl-chondrodite, hydroxyl-clinohumite and 10Å-phase. *Amer. J. Sci.* **1977**, *277*, 288–312. [CrossRef]

7. Wunder, B.; Schreyer, W. Metastability of the 10-Å phase in the system MgO-SiO₂-H₂O (MSH). What about hydrous MSH phases in subduction zones? *J. Petrol.* **1992**, *33*, 877–889. [CrossRef]

8. Comodi, P.; Cera, F.; Dubrovinsky, L.; Nazzareni, S. The high-pressure behaviour of the 10Å phase: A spectroscopic and diffractometric study up to 42GPa. *Earth Planet. Sci. Lett.* **2006**, *246*, 444–457. [CrossRef]

9. Pawley, A.R.; Welch, M.D.; Lennie, A.R.; Jones, R.L. Volume behavior of the 10-Å phase at high pressures and temperatures, with implications for H₂O content. *Amer. Mineral.* **2010**, *95*, 1671–1678. [CrossRef]

10. Bauer, J.F.; Sclar, C.B. The "10Å phase" in the system MgO-SiO₂-H₂O. *Amer. Mineral.* **1981**, *66*, 576–585.

11. Comodi, P.; Fumagalli, P.; Nazzareni, S.; Zanazzi, P.F. The 10 Å phase: Crystal structure from single-crystal X-ray data. *Amer. Mineral.* **2005**, *90*, 1012–1016. [CrossRef]

12. Momma, K.; Izumi, F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data. *J. Appl. Crystallogr.* **2011**, *44*, 1272–1276. [CrossRef]

13. Welch, M.D.; Pawley, A.R.; Ashbrook, S.E.; Mason, H.E.; Phillips, B.L. Si vacancies in the 10-Å phase. *Amer. Mineral.* **2006**, *91*, 1707–1710. [CrossRef]

14. Phillips, B.L.; Mason, H.E.; Guggenheim, S. Hydrogen bonded silanols in the 10A phase: Evidence from NMR spectroscopy. *Amer. Mineral.* **2007**, *92*, 1474–1485. [CrossRef]

15. Chen, Z.; Qi, C.; Teng, X.; Zhou, B.; Wang, C. How polar hydroxyl groups affect surface hydrophobicity on model talc surfaces. *Comm. Theor. Phys.* **2021**, *73*, 115501. [CrossRef]

16. Engelhardt, L.G.; Zeigan, D.; Jancke, H.; Wieker, W.; Hoebbel, D. ²⁹Si-NMR-Spektroskopie an Silicatlosungen. II. Zur Abhangigkeit der Struktur der Silicatanionen in wassrigen Natriumsilicatlosungen vom Na: Si-Verhaltnis. *Z. Fur Anorg. Allg. Chem.* **1975**, *418*, 17–28. [CrossRef]

17. Lippmaa, E.; Maegi, M.; Samoson, A.; Engelhardt, G.; Grimmer, A.R. Structural studies of silicates by solid-state high-resolution silicon-29 NMR. *J. Amer. Chem. Soc.* **1980**, *102*, 4889–4893. [CrossRef]

18. Cygan, R.T.; Liang, J.-J.; Kalinichev, A.G. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. *J. Phys. Chem. B* **2004**, *108*, 1255–1266. [CrossRef]

19. Cygan, R.T.; Greathouse, J.A.; Heinz, H.; Kalinichev, A.G. Molecular models and simulations of layered materials. *J. Mater. Chem.* **2009**, *19*, 2470–2481. [CrossRef]

20. Pitman, M.C.; van Duin, A.C.T. Dynamics of confined reactive water in smectite clay-zeolite composites. *J. Amer. Chem. Soc.* **2012**, *134*, 3042–3053. [CrossRef]

21. Heinz, H.; Lin, T.J.; Mishra, R.K.; Emami, F.S. Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: The INTERFACE force field. *Langmuir* **2013**, *29*, 1754–1765. [CrossRef]

22. Cygan, R.T.; Greathouse, J.A.; Kalinichev, A.G. Advances in clayff molecular simulation of layered and nanoporous materials and their aqueous interfaces. *J. Phys. Chem. C* **2021**, *125*, 17573. [CrossRef]

23. Kalinichev, A.G. Atomistic modeling of clays and related nanoporous materials with ClayFF force field. In *Computational Modeling in Clay Mineralogy; Sainz-Díaz, C.I., Ed.; Association Internationale pour l'Etude des Argiles (AIPEA): Bari, Italy, 2021; Volume 3*, pp. 17–52.

24. Underwood, T.R.; Bourg, I.C. Large-scale molecular dynamics simulation of the dehydration of a suspension of smectite clay nanoparticles. *J. Phys. Chem. C* **2020**, *124*, 3702–3714. [CrossRef]

25. Pouvreau, M.; Greathouse, J.A.; Cygan, R.T.; Kalinichev, A.G. Structure of hydrated gibbsite and brucite edge surfaces: DFT results and further development of the ClayFF classical force field with metal-O-H angle bending terms. *J. Phys. Chem. C* **2017**, *121*, 14757–14771. [CrossRef]

26. 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. [CrossRef]

27. Ho, T.A.; Greathouse, J.A.; Wang, Y.; Criscenti, L.J. Atomistic structure of mineral nano-aggregates from simulated compaction and dewatering. *Sci. Rep.* **2017**, *7*, 15286. [CrossRef] [PubMed]
28. Tararushkin, E.V.; Pisarev, V.V.; Kalinichev, A.G. Atomistic simulations of ettringite and its aqueous interfaces: Structure and properties revisited with the modified ClayFF force field. *Cem. Concr. Res.* **2022**, *156*, 106759. [CrossRef]

29. Tararushkin, E.V.; Pisarev, V.V.; Kalinichev, A.G. Equation of state, compressibility, and vibrational properties of brucite over wide pressure and temperature ranges: Atomistic computer simulations with the modified ClayFF classical force field. *Minerals* **2023**, *13*, 408. [CrossRef]

30. Fumagalli, P.; Stixrude, L. The 10Å phase at high pressure by first principles calculations and implications for the petrology of subduction zones. *Earth Planet. Sci. Lett.* **2007**, *260*, 212–226. [CrossRef]

31. Larentzos, J.P.; Greathouse, J.A.; Cygan, R.T. An ab initio and classical molecular dynamics investigation of the structural and vibrational properties of talc and pyrophyllite. *J. Phys. Chem. C* **2007**, *111*, 12752–12759. [CrossRef]

32. Molina-Montes, E.; Donadio, D.; Hernández-Laguna, A.; Sainz-Díaz, C.I.; Parrinello, M. DFT research on the dehydroxylation reaction of pyrophyllite 1. First-principle molecular dynamics simulations. *J. Phys. Chem. B* **2008**, *112*, 7051–7060. [CrossRef] [PubMed]

33. Churakov, S.V.; Dähn, R. Zinc adsorption on clays inferred from atomistic simulations and EXAFS spectroscopy. *Environ. Sci. Technol.* **2012**, *46*, 5713–5719. [CrossRef] [PubMed]

34. Hernández-Haro, N.; Ortega-Castro, J.; Pérez del Valle, C.; Muñoz-Santiburcio, D.; Sainz-Díaz, C.I.; Hernández-Laguna, A. Computational study of the elastic behavior of the 2M1 muscovite-paragonite series. *Amer. Mineral.* **2013**, *98*, 651–664. [CrossRef]

35. Escamilla-Roa, E.; Hernández-Laguna, A.; Sainz-Díaz, C.I. Theoretical study of the hydrogen bonding and infrared spectroscopy in the cis-vacant polymorph of dioctahedral 2:1 phyllosilicates. *J. Mol. Model.* **2014**, *20*, 2404. [CrossRef] [PubMed]

36. Mutisya, S.M.; Kalinichev, A.G. Carbonation reaction mechanisms of portlandite predicted from enhanced Ab Initio molecular dynamics simulations. *Minerals* **2021**, *11*, 509. [CrossRef]

37. Peng, Y.; Mookherjee, M.; Hermann, A.; Manthilake, G.; Mainprice, D. Anomalous elasticity of talc at high pressures: Implications for subduction systems. *Geosci. Front.* **2022**, *13*, 101381. [CrossRef]

38. Churakov, S.V.; Schliemann, R. Ab initio simulations of clay minerals reactivity and thermodynamics. In *Computational Modeling in Clay Mineralogy*; Sainz-Díaz, C.I., Ed.; Association Internationale pour l'Etude des Argiles (AIPEA): Bari, Italy, 2021; Volume 3, pp. 175–210.

39. Bridgeman, C.H. Ab-initio total energy study of uncharged 2:1 clays and their interaction with water. *Mol. Phys.* **1996**, *89*, 879–888. [CrossRef]

40. 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. Phy.* **2009**, *130*, 134713. [CrossRef]

41. Wang, J.; Kalinichev, A.G.; Kirkpatrick, R.J. Molecular modeling of the 10-Å phase at subduction zone conditions. *Earth Planet. Sci. Lett.* **2004**, *222*, 517–527. [CrossRef]

42. Wang, J.; Kalinichev, A.G.; Kirkpatrick, R.J. Structure and decompression melting of a novel, high-pressure nanoconfined 2-D ice. *J. Phys. Chem. B* **2005**, *109*, 14308–14313. [CrossRef]

43. Berendsen, H.; Grigera, J.; Straatsma, T. The missing term in effective pair potentials. *J. Phys. Chem.* **1987**, *91*, 6269–6271. [CrossRef]

44. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. *J. Comp. Phys.* **1995**, *117*, 1–19. [CrossRef]

45. 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. Comm.* **2022**, *271*, 108171. [CrossRef]

46. Allen, M.P.; Tildesley, D.J. *Computer Simulation of Liquids*, 2nd ed.; Oxford University Press: New York, NY, USA, 2017. [CrossRef]

47. Kühne, T.D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V.V.; Seewald, P.; Stein, F.; Guidon, M.; Holmberg, N.; Schade, R.; Krack, M.; et al. CP2K: An electronic structure and molecular dynamics software package—Quickstep: Efficient and accurate electronic structure calculations. *J. Chem. Phys.* **2020**, *152*, 194103. [CrossRef]

48. Goedecker, S.; Teter, M.; Hütter, J. Separable dual-space Gaussian pseudopotentials. *Phys. Rev. B* **1995**, *54*, 1703–1710. [CrossRef]

49. Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. *Phys. Rev. B* **1998**, *58*, 3641–3662. [CrossRef]

50. Krack, M. Pseudopotentials for H to Kr optimized for gradient-corrected exchange-correlation functionals. *Theor. Chem. Acc.* **2005**, *114*, 145–152. [CrossRef]

51. VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. *J. Chem. Phys.* **2007**, *127*, 114105. [CrossRef]

52. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868. [CrossRef]

53. 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. [CrossRef]

54. Schlegel, H.B. Optimization of equilibrium geometries and transition structures. *J. Comp. Chem.* **1982**, *3*, 214–218. [CrossRef]

55. Zanazzi, P.F.; Comodi, P.; Nazzareni, S.; Rotiroti, N.; van Smaalen, S. Behavior of 10-Å phase at low temperatures. *Phys. Chem. Miner.* **2006**, *34*, 23–29. [CrossRef]

56. Gillan, M.J.; Alfè, D.; Michaelides, A. Perspective: How good is DFT for water? *J. Chem. Phys.* **2016**, *144*, 130901. [CrossRef] [PubMed]
57. Malenkov, G.G.; Tytik, D.L.; Zheligovskaya, E.A. Hydrogen bonds in computer-simulated water. *J. Mol. Liq.* **1999**, *82*, 27–38.
[CrossRef]

58. Chanda, J.; Bandyopadhyay, S. Hydrogen bond lifetime dynamics at the interface of a surfactant monolayer. *J. Phys. Chem. B* **2006**, *110*, 23443–23449. [CrossRef]

59. Antipova, M.L.; Petrenko, V.E. Hydrogen bond lifetime for water in classic and quantum molecular dynamics. *Russ. J. Phys. Chem. A* **2013**, *87*, 1170–1174. [CrossRef]

60. Laage, D.; Stirnemann, G.; Sterpone, F.; Rey, R.; Hynes, J.T. Reorientation and allied dynamics in water and aqueous solutions. *Ann. Rev. Phys. Chem.* **2011**, *62*, 395–416. [CrossRef]

61. Yadav, S.; Chandra, A. Solvation shell of the nitrite ion in water: An ab initio molecular dynamics study. *J. Phys. Chem. B* **2020**, *124*, 7194–7204. [CrossRef]

62. Szczerba, M.; Kuligiewicz, A.; Derkowski, A.; Gionis, V.; Chryssikos, G.D.; Kalinichev, A.G. Structure and dynamics of watersmectite interfaces: Hydrogen bonding and the origin of the sharp O-DW/O-HW infrared band from molecular simulations. *Clays Clay Miner.* **2016**, *64*, 452–471. [CrossRef]

63. Blanchard, M.; Méheut, M.; Delon, L.; Poirier, M.; Micoud, P.; Le Roux, C.; Martin, F. Infrared spectroscopic study of the synthetic Mg–Ni talc series. *Phys. Chem. Miner.* **2018**, *45*, 843–854. [CrossRef]

64. Pawley, A.R.; Welch, M.D. Further complexities of the 10 Å phase revealed by infrared spectroscopy and X-ray diffraction. *Amer. Miner.* **2014**, *99*, 712–719. [CrossRef]

65. Comodi, P.; Cera, F.; Nazzareni, S.; Dubrovinsky, L. Raman spectroscopy of the 10-Å phase at simultaneously HP-HT. *Eur. J. Mineral.* **2007**, *19*, 623–629. [CrossRef]

66. Mazo, M.A.; Manevitch, L.I.; Gusarova, E.B.; Shamaev, M.Y.; Berlin, A.A.; Balabaev, N.K.; Rutledge, G.C. Molecular dynamics simulation of thermomechanical properties of montmorillonite crystal. 1. Isolated clay nanoplate. *J. Phys. Chem. B* **2008**, *112*, 2964–2969. [CrossRef] [PubMed]

67. Mazo, M.A.; Manevitch, L.I.; Gusarova, E.B.; Berlin, A.A.; Balabaev, N.K.; Rutledge, G.C. Molecular dynamics simulation of thermomechanical properties of montmorillonite crystal. II. Hydrated montmorillonite crystal. *J. Phys. Chem. B* **2008**, *112*, 17056–17062. [CrossRef]

68. Hill, R. The elastic behaviour of a crystalline aggregate. *Proc. Phys. Soc. Sect. A* **1952**, *65*, 349–354. [CrossRef]

69. Mainprice, D.; Page, Y.L.; Rodgers, J.; Jouanna, P. Ab initio elastic properties of talc from 0 to 12 GPa: Interpretation of seismic velocities at mantle pressures and prediction of auxetic behaviour at low pressure. *Earth Planet. Sci. Lett.* **2008**, *274*, 327–338. [CrossRef]

70. Kostenetskiy, P.S.; Chulkevich, R.A.; Kozyrev, V.I. HPC resources of the Higher School of Economics. *J. Phys. Conf. Ser.* **2021**, *1740*, 012050. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.