Insight into the Structure of TMA-Hectorite: A Theoretical Approach

: An ab initio density functional theory method (DFT) with D3 dispersion corrections (DFT-D3) was employed to study the 64 possible models of the structure of hectorite intercalated with tetramethylammonium (TMA) cations with the aim to ﬁnd the additional information about the structure for the energetically most suitable mutual arrangement of the TMA cations. The complex analysis of total energies showed small differences among the structural models with the lowest (L model) and the highest (H model) total energy (~50 kJ/mol). The analysis of the calculated vibrational spectrum of the L model of the TMA-H structure was conducted in detail.


Introduction
Natural or synthetic smectites, like hectorites, are important, widely abundant and low-cost materials with unique physicochemical properties, such as swelling, intercalation and ion-exchange properties. Synthetic hectorites are especially applicable, e.g., in catalytic processes [1] or matrices in polymer nanocomposites, because of their controlled pore size distribution and composition in comparison to natural minerals [2]. The characteristics mentioned above also predetermine the hectorites for application in medicine and pharmacology, e.g., for drug delivery purposes [3,4] or as a potential sorbent of waste compounds [5]. In spite of the wide use of the smectites (also hectorites) in the industry, the structural characterization of clays and afterwards organoclays is very difficult because of a poor crystallinity of these materials. The information obtained from the powder X-ray measurements are very often insufficient. The positions of hydrogen atoms are missing, hence the information about hydrogen bonds of the intercalated molecules is lacking, except for quite rare experiments based on neutron powder diffraction [6][7][8]. To investigate the clay surface and to characterize the interactions between organic compounds and clay minerals, some experimental methods, such as Fourier transform infrared spectroscopy (FTIR) [9], Ultraviolet-visible (UV-vis) spectroscopy [10], contact angle measurement [11], atomic force microscopy [12], NMR [13] and He-ion imaging experiments [14] have been employed.
The computational methods help to eliminate the deficiency of suitable information about organoclays structures, to complete the experimental data and to find more and more applications in the fields of organoclays [15][16][17][18][19][20][21][22][23][24][25] and mineralogy [26][27][28][29][30]. In addition to the refined atomic positions, e.g., hydrogen atoms, which enable the study of hydrogen bonds, it is also possible to analyze the vibrational modes of the individual functional groups in these structures and to help with analyzing the overlapped vibrational bands in the experimentally measured spectrum [31][32][33][34]. When using computational methods, the proposing of a suitable representative model of the organoclay structure is an important factor. In addition to the charge of the intercalated organic cations in the interlayer space, the substitutions in the clay structure also have to be taken into consideration during the construction of the model [35]. The determination of the atomic positions of the intercalated molecules/cations might be problematic as well [36], because the further experimental refinement of the atomic positions of the intercalated species is not possible. This is also the

Computational Details
DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) program [41,42]. The exchange-correlation energy was expressed in the framework of the generalized gradient approximation (GGA) using the DFT functional proposed by Perdew, Burke and Ernzerhof (PBE) [43]. The Kohn-Sham equations were solved using a variational method in a plane-wave (PW) basis set with an energy cut-off of 500 eV. The electron-ion interactions were described using the projector-augmented-wave (PAW) method [44,45]. Brillouin-zone sampling was restricted to the gamma point. The atomic positions were relaxed in the restricted experimentally obtained unit cell parameters. The relaxation criteria were 10 −5 eV/atom for the total energy change and 0.005 eV/Å for the maximum force acting on any atom.
Normal modes of vibrations were calculated within the fixed optimized cell using a finite difference method and harmonic approximation. The Hessian was constructed from single point energy calculations on the 6N structures generated from the optimized structures by displacing each of the N th atoms in the cell in a positive and negative direction along the three Cartesian directions x, y, and z [46]. DFT calculations were performed involving dispersion corrections using a D3 scheme, as recommended for the PBE functional [47].

Computational Models
Computational models were based on the structural data suggesting two possible orientations of TMA cations with similarly low occupancies (0.09 and 0.08) [38]. It is evident that with such a low value, it is not possible to keep a reasonable size of the computational cell and mimic the refined occupancies at the same time. The proposed models are thus based on the assumption that in a computational (a, b, c) cell there are, along with the 2:1 layer, two TMA cations, whose central nitrogen atoms could only appear in any of the positions constrained in the original space group (C2/m) to the 8j Wyckoff position, WP (Table 1). In the P1 space group, they are, however, free to translate and rotate, though all of that movement is restricted by the potential field generated by the surrounding atoms. In order to sample the space of mutual orientations, a simple believably sufficient strategy was suggested. The strategy of positioning of the cations labeled here as A and B could best be explained the by using an 8 × 8 table made of four 4 × 4 quadrants. The design of such a table is provided in Table 2 where the "instructions" on how to generate the individual configurations are available. For instance, if only one type (A or B) of TMA cations is to appear in the upper left quadrant of the 8x8 configurations table, one takes X = Y = A and combines it with the sequential number of the symmetry operation of the general 8j WP. The configurations based on the B cation are in the lower right part of the configurations table, i.e., X = Y = B. Quite naturally, the lower left corner is generated by taking X = B and Y = A, while in the upper right one the order is reversed, i.e., X = A and Y = B. The auxiliary labels for respective A (a 1 -a 8 ) and B (b 1 -b 8 ) TMA-cations were added for better clarity, i.e., which X, Y positions were used for A and/or B types of TMA cation in the respective model (Table 1). In such a way, all 64 computational models were formed. In each case, the  Finally, three optimized models were selected for the next considerations: L model with the lowest total energy (Figure 1), H model with the highest energy and R model as the randomly chosen energy, representing the group of configurations with close energies higher than in the L model and lower than in the H model.  Finally, three optimized models were selected for the next considerations: L model with the lowest total energy (Figure 1), H model with the highest energy and R model as the randomly chosen energy, representing the group of configurations with close energies higher than in the L model and lower than in the H model. , (x + 0.5, y + 0.5, z) Wyckoff position of the second TMA cation (a5 for A type) as an example (Si-tetrahedra in yellow, Mg-octahedra in cyan and Li-octahedra in green. N atoms in blue, C-atoms in light grey, O atoms in red, H atoms in white and F atoms in light green). Table 1. X, Y Wyckoff position, 8j, for atoms in A and/or B types of TMA cation, respectively.

Structure and Total Energy Analysis
Atomic coordinates in 64 models were fully optimized and graphical distribution of the resulting total energies is presented in Figure 2. There are only very small differences in energy among the majority of the models. The difference between the highest and Minerals 2021, 11, 505 4 of 11 lowest energy was~0.5 eV (~50 kJ/mol). Approximately 70% of configurations (43 models) depicted in the histogram (Figure 3, down) fell into the region with very close energies and differentiation between configurations is not easy, if even possible. This fact was also confirmed by cumulative counts' graph ( Figure 3, up) where the percentage of the counts for respective bar (including all bars with smaller values of the response variable) is presented. The respective percentage of cumulative counts for the bar with 43 models was 85% (55 models). Three models were selected for the next considerations as mentioned in the previous paragraph: L model with the lowest total energy, −451.248 eV, WP of N atom of TMA (1) cation is 4X (a 4 for A type) and for TMA (2)  Atomic coordinates in 64 models were fully optimized and graphical distribu the resulting total energies is presented in Figure 2. There are only very small diffe in energy among the majority of the models. The difference between the highest an est energy was ~0.5 eV (~50 kJ/mol). Approximately 70% of configurations (43 m depicted in the histogram (Figure 3, down) fell into the region with very close en and differentiation between configurations is not easy, if even possible. This fact w confirmed by cumulative counts' graph ( Figure 3, up) where the percentage of the for respective bar (including all bars with smaller values of the response variable) sented. The respective percentage of cumulative counts for the bar with 43 mode 85% (55 models). Three models were selected for the next considerations as mentio the previous paragraph: L model with the lowest total energy, −451.248 eV, WP of N of TMA (1) cation is 4X (a4 for A type) and for TMA (2)     The analysis of the optimized structural parameters showed that the bond distances within the 2:1 layer were very close to experimental data reported by Seidl and Breu (2005) [38] for all selected models. The optimized bond distances for the L model, as an example, compared to experimental data differed in the hundredths of Å (Table 4).  The analysis of the optimized structural parameters showed that the bond distances within the 2:1 layer were very close to experimental data reported by Seidl and Breu (2005) [38] for all selected models. The optimized bond distances for the L model, as an example, compared to experimental data differed in the hundredths of Å (Table 4).

Hydrogen Bonds
Calculated C···O (Donnor···Acceptor) distances for the TMA cations for selected models and the distances published by Seidl and Breu (2005) [38] are given in Table 5. These C···O contact distances belong to the weak hydrogen bonds [48,49]. A comparison of the hydrogen bonds among the models showed that their character and strength did not differ much, as should be expected. The differences were mainly in the mutual orientation of the both TMA cations, as the numbers and strength of hydrogen bonds slightly decreased with the total energy of the respective models (Figure 4a-c, Tables 3 and 5). Some hydrogen bonds were two centered. Total energy as a measure of stability of the models correlated well with the values of the C-H···O (Donnor-Hydrogen···Acceptor) angles in the hydrogen bonds. In general, the sharper the angle was, the weaker the hydrogen bond became. The strength of the bonds decreased from the L model towards the energetically less stable H model (Table 3). the both TMA cations, as the numbers and strength of hydrogen bonds slightly decreased with the total energy of the respective models (Figure 4a-c, Tables 3 and 5). Some hydrogen bonds were two centered. Total energy as a measure of stability of the models correlated well with the values of the C-H···O (Donnor-Hydrogen···Acceptor) angles in the hydrogen bonds. In general, the sharper the angle was, the weaker the hydrogen bond became. The strength of the bonds decreased from the L model towards the energetically less stable H model (Table 3).

Powder Diffraction Patterns
Theoretical powder diffraction patterns were also calculated for the L and H models in order to find possible differences reflecting the changes in the mutual orientations of TMA cations in the models. Again, one should not expect that different mutual positions of TMA cations produce very different energies of the structures. Considering that the majority of the TMA models give very similar energies, the changes in the powder diffraction patterns should be very small as well. This hypothesis was confirmed, as can be seen in Figure 5. Comparison of both calculated powder diffraction patterns for the L and H model structures showed that the differences were almost negligible and the distinguishing of these two structures is impossible.

Powder Diffraction Patterns
Theoretical powder diffraction patterns were also calculated for the L and H models in order to find possible differences reflecting the changes in the mutual orientations of TMA cations in the models. Again, one should not expect that different mutual positions of TMA cations produce very different energies of the structures. Considering that the majority of the TMA models give very similar energies, the changes in the powder diffraction patterns should be very small as well. This hypothesis was confirmed, as can be seen in Figure  5. Comparison of both calculated powder diffraction patterns for the L and H model structures showed that the differences were almost negligible and the distinguishing of these two structures is impossible.

TMA Cations Anchoring
To compare the optimized atomic positions of the TMA cations in the L model with experimental data, the optimized atomic positions of relevant atoms were transformed to

TMA Cations Anchoring
To compare the optimized atomic positions of the TMA cations in the L model with experimental data, the optimized atomic positions of relevant atoms were transformed to the initial x, y, z WP atomic position (Table 1) published by Seidl and Breu (2005) [38]. Comparison of the optimized atomic positions with experimental data showed a good agreement. The differences among the experimental and optimized atomic positions of respective atoms were very small for both the hectorite structure (Table 6) and TMA cations ( Table 7). The optimized atomic positions of the hydrogen atoms are presented in Table 8 to complete the structural information on the TMA cations in the interlayer space of hectorite. Table 6. Optimized atomic positions in the energetically most stable L model of the TMA-H structure in comparison to the experimental structure [38].

Atoms
x  Further, the position of TMA cations in the interlayer space is similar to the one in the refined model in the work of Seidl and Breu (2005) [38]. That means that TMA cations are positioned in the middle of the interlayer space. Two of the methyl groups lie in the ab plane, whereas the remaining two methyl groups point up and down ( Figure 6a) and they are centered in the ditrigonal hole ( Figure 6b).  The mutual configuration of the two TMA cations in the optimized TMA-H structure (L model) is presented in Figure 7. Neighboring cations are rotated in order to avoid unrealistically close contacts. N···Obasal distances vary in the range of 3.84-4.71 Å. Both cations are anchored into the silicate tetrahedral sheets on both sides and in this way cross-link the interlayer space through hydrogen bonds discussed in the previous paragraph ( Figure  4a).

Vibrational Spectra
To complete the information about the TMA-H structure, the vibrational spectrum for the L model was calculated (Figure 8). The respective calculated unscaled vibrational bands can be formally divided into two regions: the first is 3250-2750 cm −1 and the second

Vibrational Spectra
To complete the information about the TMA-H structure, the vibrational spectrum for the L model was calculated (Figure 8). The respective calculated unscaled vibrational bands can be formally divided into two regions: the first is 3250-2750 cm −1 and the second is below 1500 cm −1 .

Conclusions
An ab initio DFT-D3 method was employed for examination of possible mutual orientations (64 models) of TMA cations in the hectorite interlayer space (TMA-H) with the aim to supply new information about the structure of TMA-H organoclays, such as a final refinement of hydrogen atoms positions. From systematical study of 64 models of the TMA-H modified clay, only very small differences in total energy among the majority of the models were found (~50 kJ/mol as the maximum). Powder diffraction patterns of the models with the highest (H) and lowest energy (L) were not distinguishable. The strength of the weak C···O hydrogen bonds presented in the interlayer space of the hectorite slightly decreases with the total energy of the studied systems. Calculated vibrational modes for the energetically most stable L model of TMA-H organoclay were discussed in detail.  Three main bands are present in the first region and one of them has a shoulder. The detailed analysis of the individual calculated vibrational modes showed that the highest energy band in the calculated spectrum of the L model corresponds with the C-H asymmetric stretching vibrations, with a maximum at 3136 cm −1 . The C-H symmetric stretching vibrations are present in the rest of two bands, with the maximum energy at 2995 cm −1 (shoulder at 3034 cm −1 ) and at 2927 cm −1 .
In the second region, the highest energy band is an association of the two bands with a maximum at 1450 cm −1 for C-H asymmetric bending vibration and at 1410 cm −1 for common vibrations of C-N (asymmetric stretching) and C-H (asymmetric bending) groups. The C-N symmetric stretching appears at 1287 cm −1 and C-H symmetric bending at 1240 cm

Conclusions
An ab initio DFT-D3 method was employed for examination of possible mutual orientations (64 models) of TMA cations in the hectorite interlayer space (TMA-H) with the aim to supply new information about the structure of TMA-H organoclays, such as a final refinement of hydrogen atoms positions. From systematical study of 64 models of the TMA-H modified clay, only very small differences in total energy among the majority of the models were found (~50 kJ/mol as the maximum). Powder diffraction patterns of the models with the highest (H) and lowest energy (L) were not distinguishable. The strength of the weak C···O hydrogen bonds presented in the interlayer space of the hectorite slightly decreases with the total energy of the studied systems. Calculated vibrational modes for the energetically most stable L model of TMA-H organoclay were discussed in detail.
Data Availability Statement: Important structural data are presented directly in the manuscript.