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

: 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 3 Si 4 O 10 (OH) 2 · x H 2 O, 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 2 O molecules within the mineral are not yet sufﬁciently 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 2 O 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 modiﬁcation of the ClayFF force ﬁeld 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 2 O molecules.


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][2][3][4][5]. Therefore, TAP structural, thermodynamic, and mechanical properties are important, especially at high pressures and temperatures [6][7][8][9]. TAP is assumed to have the composition of Mg 3 Si 4 O 10 (OH) 2 ·xH 2 O 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 2 O 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 2 O 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 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 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. . Green octahedra-MgO, blue tetrahedra-SiO, red spheres-O, white spheres-H. VESTA software [12] was used to visualize the atomistic model. 29 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 2 -type Si site (see Figure S1 of the Supplementary Materials) contains silanol groups that donate hydrogen bonds (H-bonds) to the interlayer H2O molecules. These results suggest that the formation of TAP can involve a defect mechanism that allows favorable H-bonding interaction between the interlayer H2O molecules and the normally hydrophobic siloxane oxygens of the talc-like layer [14,15]. The proportion of the Q 2 -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 n , where the superscript shows the number of other silicon-oxygen tetrahedra attached to the silicon tetrahedron under study [16,17]).  [12] was used to visualize the atomistic model. 29 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 2 -type Si site (see Figure S1 of the Supplementary Materials) contains silanol groups that donate hydrogen bonds (H-bonds) to the interlayer H 2 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 2 O molecules and the normally hydrophobic siloxane oxygens of the talc-like layer [14,15]. The proportion of the Q 2 -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 n , 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 atomistic modeling can be used as a powerful tool to clarify this picture. Classical MD simulations 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][19][20][21][22][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][23][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][28][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][31][32][33][34][35][36][37][38][39][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]). hydrophobic siloxane oxygens of the talc-like layer [14,15]. The proportion of the Q 2 -Si sites inferred from the NMR data is estimated to be around 10% [14], correspondin one silanol for each pair of six-membered siloxane rings of the talc tetrahedral sheet. ( we are using the common notation for silicon-oxygen tetrahedra, Q n , where the su script shows the number of other silicon-oxygen tetrahedra attached to the silicon t hedron under study [16,17]).

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.
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.

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 2 -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.

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) 2 ), gibbsite (Al(OH) 3 ), and kaolinite (Al 2 Si 2 O 5 (OH) 4 ) using a simple harmonic functional form [22,25,26]: 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 kcal·mol −1 ·rad −2 , θ 0,SiOH = 100 • , and k SiOH = 15 kcal·mol −1 ·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.

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][49][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 9 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 −6 a.u. All AIMD Minerals 2023, 13, 1018 5 of 17 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.

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 1. Crystallographic unit cell parameters of the TAP models: Simulated results and experimental data. 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 aand 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 2 O molecules in the structure [25,56].

Atom Positions and Interatomic Distances
All calculated average distances between metal atoms (Mg, Si) and various oxygen atoms in the structure (O a -apical oxygen atom; O b -basal oxygen; O h -oxygen atom of hydroxyl; O w -oxygen atom of H 2 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. 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 s and O bs (Si and basal O belonging to silanol groups) was slightly larger than the Si-O b distance in all simulations. The same result was also obtained for the Si-O a pair in the AIMD calculations ( 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.

Radial Distribution Functions
Radial distribution functions (RDFs) are especially important for understanding the structural arrangement and ordering of the H 2 O molecules in the TAP interlayers. They were calculated for different H-O pairs existing in our models. The first H w -O w 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 2 O molecules compared to the experimental data ( Table 2). water molecules. This effect was not observed with the ClayFF-orig force field.
The position of the first Hw-Ob 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 Hw-Ob 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 Hh-Ow 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 H2O molecules in the six-membered tetrahedral rings [11], which were free from short-range interaction due to the structural defects. The first peak for the ideal model was slightly shifted to the left in the CMD simulations, 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 hs ) and O w 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 2 O molecules approaching the tetrahedral layer. Consequently, 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 w -O b 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 w -O b 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 h -O w 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 2 O molecules in the six-membered tetrahedral rings [11], which were free from short-range interaction due to the structural defects.  Fairly good agreement with the ClayFF-MOH and AIMD simulation results was also 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 hs ···O w 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 hd -O bs RDFs ( Figure 5; see also Figure 2). The average distance between H hd and O bs 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 bs atoms for H-bond formation with one of the H hd 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 hd and O bs was most likely.  This suggestion was also confirmed by the Hhs-Ow 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 Hhs···Ow 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 Hhd-Obs RDFs ( Figure 5; see also Figure 2). The average distance between Hhd and Obs 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 Obs atoms for H-bond formation with one of the Hhd 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 Hhd and Obs was most likely.

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).
For the ideal TAP model, there were no differences between the two versions of ClayFF and AIMD results. In all cases, the H2O molecules were located above the structural hydroxyls of each six-membered siloxane ring [41] of the talc tetrahedral layers (see Figure 7a). Also, the H2O molecules were coordinated by the Ob atoms of the tetrahedral layer [41]. According to previous DFT results [30], the Hw atom can form multi-furcated H-bonds with Ob atoms.
Water molecules located far from the structural defects had similar behaviors to the H2O molecules in the ideal model for both ClayFF versions ( Figure 6). However, the H2O 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 H2O molecule arrangement near the defects, which was due to their acceptance of strong H-bonds from the silanol groups of the defects (Figure 7b).

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).

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]: ROdOa ≤ 3.5 Å, ROaHd ≤ 2.45 Å, and φHdOdOa ≤ 30°, where d is a donor and a is an acceptor. The average lifetime of H-bonds (τHB) was estimated by integrating the so-called continuous time functions of H-bonds [58,59]. The average number of H-bonds (nHB) per donor was also calculated (Table 3;   For the ideal TAP model, there were no differences between the two versions of ClayFF and AIMD results. In all cases, the H 2 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 2 O molecules were coordinated by the O b atoms of the tetrahedral layer [41]. According to previous DFT results [30], the H w atom can form multi-furcated H-bonds with O b atoms.
Water molecules located far from the structural defects had similar behaviors to the H 2 O molecules in the ideal model for both ClayFF versions ( Figure 6). However, the H 2 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 2 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).

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 OdOa ≤ 3.5 Å, R OaHd ≤ 2.45 Å, and ϕ HdOdOa ≤ 30 • , where d is a donor and a is an acceptor. The average lifetime of Hbonds (τ HB ) was estimated by integrating the so-called continuous time functions of H-bonds [58,59]. The average number of H-bonds (n HB ) per donor was also calculated (   The short lifetimes of H-bonds indicated that the H 2 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 w -H w bond in the H 2 O molecules ( Figure 8). Thus, the interaction between structural hydroxyls and H 2 O molecules played a very important role in the behavior of the interlayer H 2 O molecules in the ideal model.  In the new TAP model with defects, several new donor-acceptor H-bonding pairs were possible in addition to the main Hw···Ob and Hh···Ow H-bonds. The most important was the Hhs···Ow pair. It had the longest lifetime, especially as observed in the AIMD simulations ( Table 3 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 donoracceptor pair. This was due to the possibility of frequent and completely unrestricted reorientation of the silanol hydroxyl groups. The average number of Hhs···Ow 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 Hw···Obs, which was due to H-bond donation by the H2O 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

Defect Model
The observed behavior of the H 2 O molecules far from the TAP structural defects was the same as in the ideal model. The H w ···O b H-bonding lifetime did not change significantly (Table 4). However, the number of H w ···O b H-bonds was smaller in all simulations. On the contrary, an increase in the average number of H h ···O w H-bonds was observed for the ClayFF-MOH and AIMD simulations. The lifetime of H h ···O w 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 2 O molecules. The orientational relaxation time was the 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 2 O 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 2 O molecules played even a greater role in the defect TAP model.

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 2 O 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]. The stretching modes of the H2O 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 4). A more accurate and elaborate H2O intramolecular potential would be necessary to better reproduce the positions and widths of these peaks (e.g., [62]). The Ow-Hw 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.
The AIMD peak positions reflecting the stretching mode of structural hydroxyls (Hh atoms) were blue-shifted compared to the experimental data and the results of the CMD The stretching modes of the H 2 O 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 2 O 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. The AIMD peak positions reflecting the stretching mode of structural hydroxyls (H h 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 −1 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 hs atoms) (Figure 9). The peak position obtained from the AIMD simulations was located at 3194 cm −1 , while the CMD results showed peaks located at 3577 and 3431 cm −1 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 2 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 redshifted compared to the high-pressure experimental data of Pawley and Welch [64], in which these vibrational modes were detected at 3587 cm −1 . 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 h 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 2 O were located at ca. 310−320 cm −1 in the CMD simulations and at ca. 375−385 cm −1 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.

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. 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.

Conclusions
Six atomistic models of the 1M polytype of the 10-Å phase (TAP), Mg 3 Si 4 O 10 (OH) 2 ·xH 2 O, were constructed and quantitatively studied using classical and ab initio molecular dynam-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 2 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 2 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 3 -and Q 2 -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 2 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.