# Structure and properties of thaumasite and its aqueous interfaces revealed by molecular dynamics simulations

Evgeny V. Tararushkin, Grigory S Smirnov, Nikolay D. Kondratyuk, Andrey

G. Kalinichev

► To cite this version:

Evgeny V. Tararushkin, Grigory S Smirnov, Nikolay D. Kondratyuk, Andrey G. Kalinichev. Structure and properties of thaumasite and its aqueous interfaces revealed by molecular dynamics simulations. Cement and Concrete Research, 2025, 197, pp.107944. 10.1016/j.cemconres.2025.107944. hal-05153734

HAL Id: hal-05153734

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

Submitted on 14 Jul 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>

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

# Structure and properties of thaumasite and its aqueous interfaces revealed by molecular dynamics simulations

E.V. Tararushkin<sup>a, ID,*</sup>, G.S. Smirnov<sup>a, ID</sup>, N.D. Kondratyuk<sup>b, a ID</sup>, A.G. Kalinichev<sup>c, ID,*</sup>

## ARTICLE INFO

MSC:

0000

1111

Keywords:

Thaumasite

Ettringite

Hydrogen bonds

Interfaces

Interfacial energy

Epitaxial growth

Molecular dynamics

DFT

PIMD

## ABSTRACT

Thaumasite, along with ettringite, is responsible for sulfate attack in concrete based on Portland cement. We use classical molecular dynamics (MD) computer simulation technique together with the ClayFF-MOH force field to quantitatively investigate the structural, vibrational, mechanical and thermodynamic properties of thaumasite. Additionally, path integral MD technique is applied for a more detailed examination of the structure and dynamics of intracrystalline hydrogen bonds in this material. The results of classical MD simulations are in good agreement with available experimental data and DFT calculations. Mineral-water interfaces for the (1 0 0) and (0 0 1) surfaces of thaumasite and ettringite are also simulated and compared. The solid-water interfacial energies for (1 0 0) and (0 0 1) thaumasite and ettringite surfaces provide new insights on the atomistic scale solubility mechanisms of these minerals. The simulations also confirm the possibility of thaumasite epitaxial growth on the surface of ettringite.

## 1. Introduction

Thaumasite is a hydrated calcium silanol-carbonate-sulfate natural mineral occurring in sulfide and carbonate metamorphic rocks [1–3]. Chemically represented as Ca₃[Si(OH)₆](CO₃)(SO₄)·12H₂O [2,4,5], thaumasite exhibits a columnar crystal structure similar to that of ettringite (Ca₆[Al(OH)₆]₂(SO₄)₃·nH₂O), which is comprised of stacked silicon and calcium hydroxide polyhedra [Ca₃Si(OH)₆(H₂O)₁₂]⁴⁺. The hydroxide columns of thaumasite running parallel to the c lattice vector [2,6], are interconnected through a network of hydrogen bonds (H-bonds), facilitated by H₂O molecules, SO₄²⁻ and CO₃²⁻ anions that occupy the intercolumnar spaces within the crystal (see Fig. S1 in the Electronic Supplementary Information). In the hydroxide columns, silicon is octahedrally coordinated by hydroxyls and this makes the structure of thaumasite extremely unusual from the crystallochemical perspective [7].

Under normal conditions, the crystal chemistry of silicon–oxygen compounds is dominated by the tetrahedral 4-fold coordination of silicon atoms, while the octahedral 6-fold coordination is only found in minerals under extreme conditions of temperature and pressure, such as existing in the Earth’s mantle [8,9]. Silica atoms in 5-fold and 6-

fold coordination can also be occasionally found in silica glass [10,11], or as an intermediate species in the process of aqueous dissolution of phyllosilicates [12]. However, thaumasite represents a rare example of a crystal in which silicon is octahedrally coordinated under ambient temperature and pressure away from these harsh conditions.

In addition to its importance in mineralogy and crystallography, thaumasite is of considerable interest within the domain of cement and concrete chemistry [13]. The thaumasite form of the sulfate attack (TSA) is a phenomenon that can cause damage to concrete at low temperature (0–5 °C in particular) due to the expansion of the cement matrix [4,5,14–17]. TSA is commonly observed in building constructions such as foundations and tunnels, where concrete is exposed to soil sulfate-containing solutions in conjunction with the carbonation process [18].

The formation of thaumasite (TF) is assumed to be preceded by decomposition of the calcium silicate hydrate (C–S–H) phase in the presence of dissolved sulfates and carbonates in the pore solution of cementitious materials in excess humidity [4,15,18–20]. However, the molecular mechanisms of TF are still unclear. Three main TF routes are considered [4,14,18]. One is the direct route from C–S–H reacting with dissolved CO₃²⁻, SO₄²⁻, Ca²⁺ ions, and excess water [14]. The second

<footnote><sup>a</sup> International Laboratory for Supercomputer Atomistic Modeling and Multi-Scale Analysis, HSE University, Moscow, 101000, Russia</footnote>

<footnote><sup>b</sup> Center for Computational Physics, Moscow Institute of Physics and Technology, Dolgoprudny, 141701, Russia</footnote>

<footnote><sup>c</sup> Laboratoire SUBATECH (UMR 6457—IMT Atlantique, Nantes Université, CNRS/IN2P3), Nantes, 44307, France</footnote>

<footnote>* Corresponding authors.</footnote>

<footnote>E-mail addresses: evgeny.tararushkin@yandex.ru (E.V. Tararushkin), kalinich@subatech.in2p3.fr (A.G. Kalinichev).</footnote>
pathway involves the formation of thaumasite in a solid solution with ettringite (known as the woodfordite route), a compound often found in cement-based materials [14]. The third path involves the epitaxial growth of thaumasite on the surface of ettringite [4].

Until recently, various experimental techniques were the main tools in studying the properties and behavior of cementitious materials. However, experimental possibilities are often limited, especially for quantitative understanding of the molecular-scale structure and processes. However, over the past two decades, computational atomistic modeling approaches have already established themselves as reliable complementary tools to quantitatively investigate the structure and properties of cementitious materials at the nano-scale. Classical and *ab initio* molecular dynamics simulations are successfully used to study cementitious materials and their aqueous interfaces [21,22]. In particular, structural, mechanical, vibrational, and electronic properties of thaumasite have been investigated earlier by molecular mechanics and molecular dynamics simulations [1,23–26]. The important role of the intracrystalline H-bonds was also analyzed [23]. However, similar investigations of the thaumasite interfaces with aqueous solutions are still lacking, even though such an analysis could greatly improve our nano-scale understanding of the TF mechanisms.

Here we are using classical molecular dynamics simulations (MD) to study the crystal properties of thaumasite, the dynamics and structural role of its intracrystalline H-bonds, and the structure of aqueous solutions at its various crystal faces. Several empirical force field models are often used in classical MD simulations of cementitious materials, and each model usually demonstrates the most reliable results in investigations of the properties and phenomena for which they have been specifically validated to reproduce the experimentally available data [22]. Although there is no single “best” force field for all possible applications and phenomena, the principal limitation of any empirical force field approach is that it cannot be used to model chemically reactive systems, which is often important in cement chemistry.

However, within its known limitations, the ClayFF force field [27] used here to describe all interatomic interactions has already proven to be a simple and robust tool for classical MD simulations of many layered, nanoporous, and nanostructured materials, including clay and cement [22,28]. Its most recent modification, ClayFF-MOH, explicitly takes into account Metal–O–H (M–O–H) angular bending motion of the hydroxyl groups in the mineral structure [29–32], leading to a much better description of the hydroxyl behavior in the crystal structures, at the edges of mineral nanoparticles, and on irregular surfaces [33, 34]. Recent atomistic simulations of ettringite and its aqueous interfaces [32] have already clearly demonstrated the importance of the M–O–H terms for accurate modeling of this mineral, which is structurally similar to thaumasite.

In addition to classical MD modeling, we also explore electronic and nuclear quantum effects. Electronic quantum effects are treated using the density functional theory (DFT) approach [35]. Nuclear quantum effects (NQEs), such as tunneling and zero-point energy, play a significant role in systems containing light atoms [36,37]. NQEs cannot be reproduced neither by classical nor by quantum MD simulations due to the classical treatment of the nuclear subsystem in both approaches. Path Integral Molecular Dynamics (PIMD) is one of the simulation techniques that incorporates the quantum nature of the nuclei [36,37]. PIMD is particularly effective for simulating water and hydrous materials (e.g., clay or cement) that contain hydrogen or lithium [38–45]. This recent research indicates that NQEs can lead to a weakening of H-bonds network due to the hydrogen delocalization. Such effects are especially important at high pressures, as has been shown for brucite and portlandite [46,47].

Finally, we investigate the structural characteristics of water on the various thaumasite surfaces using the ClayFF potential. In addition, we calculate the solid–liquid interfacial energy. All of these properties are compared with similar aqueous interfaces of ettringite.

## 2. Models and methods

### 2.1. Structural models

The structural model of thaumasite was based on the single-crystal X-ray diffraction data, which provides the following unit cell parameters: $a = 11.0538(6)$ Å, $c = 10.4111(8)$ Å, $V = 1101.67(10)$ Å³ with the $P\bar{6}_3$ space group [5] (see Fig. S1 in the Electronic Supplementary Information). The columnar crystal structure of thaumasite is similar to that of ettringite, but with approximately 1/2 of the unit cell in the c-dimension due to the higher crystal symmetry of thaumasite [48].

The simulation supercell of thaumasite containing 3294 atoms (3 × 3 × 3 unit cells) was used in the classical MD simulations. For DFT and quantum molecular dynamics (QMD) calculations, smaller supercells with only 488 atoms (2 × 2 × 1 unit cells) were used. Periodic boundary conditions were applied in all three dimensions to produce continuous models of crystals.

Similarly to ettringite [32], we studied the thaumasite surfaces in contact with water using classical MD simulations with the ClayFF-MOH/harmonic model. For these simulations, ettringite and thaumasite supercells were cleaved parallel to the (1 0 0) and (0 0 1) crystallographic planes. The positive charge of the cleaved model surfaces is formed due to the hydroxyls and the hydrogen atoms of the interfacial H₂O molecules. However, presence of the surface anions results in some regions having a negative charge. Therefore, all modeled surfaces were assumed to have a heterogeneous charge density distribution, but the actual values are still unknown. Nevertheless, it was shown that the surface charge of ettringite depends on the concentrations of positively and negatively charged species in the interfacial aqueous solution [49,50]. We assume a similar effect for thaumasite surfaces.

All cleaved surfaces were brought into contact with an approximately 80 Å thick layer of water molecules, as shown in Fig. S2 and S3. The dimensions of the resulting simulation supercells were approximately 70 × 70 × 130 Å³ for the (1 0 0) interface and 65 × 57 × 147 Å³ for the (0 0 1) interface.

### 2.2. Force field parameters and details of classical MD simulations

The standard ClayFF parametrization does not include octahedrally coordinated Si atoms [29]. Therefore, we modeled these atoms by assigning them the same Lennard-Jones (L-J) parameters as for octahedral Al ($\sigma = 4.2712$ Å, $\epsilon = 5.5638832 \times 10^{-6}$ kJ/mol), while assigning the same atomic charge as for tetrahedrally coordinated silicon $q = +2.1e$. It will be shown below that this choice of octahedral silicon parameters is quite reasonably justified. To bring this thaumasite crystal structural model to neutrality, the charge of –0.7883317e was assigned to the oxygen atoms of the structural hydroxyls. This value is within the range of typical partial atomic charges for oxygen atoms in various structural environments accepted by the ClayFF force field [29,51].

The non-bonded and bonded force field parameters for the $SO_4^{2-}$ anions were the same as in our previous study of ettringite [32], while the parameters for $CO_3^{2-}$ anions were taken from Ref. [52]. H₂O molecules were described by the extended simple point charge model (SPC/E) with flexible intramolecular O–H bonds [29,53].

For the O–H bonds stretching in H₂O molecules and in structural hydroxyls, we compare the harmonic and Morse representations. Although the harmonic potential was used in the original ClayFF parametrization [27], the more recent Morse modification allows for a more accurate description of the O–H stretching modes [29,54].

In total, we have tested three ClayFF modifications with the following notation: ClayFF-MOH/Morse, ClayFF-MOH/harmonic, and ClayFF/Morse. However, the last combination, based on the original parametrization [27] that did not explicitly take into account the M–O–H bending terms in the structure, results in an unstable thaumasite crystal and will not be discussed further. Lennard-Jones parameters
for the unlike atom types were calculated using the Lorentz–Berthelot combining rules [29].

LAMMPS software package [55] was used to perform all classi-
cal MD simulations. To speed up the calculations, we utilized the
GPU accelerators [56,57]. To account for the long-range electrostatic
interactions, the particle–particle particle–mesh method (PPPM) was
used [58] with the cutoff radius of 12.5 Å. The velocity Verlet algo-
rithm [59] was used to numerically integrate the atomic equations of
motion with the 1 fs timestep. Initially, the structures were equilibrated
in the *NPT* ensemble for 1 ns at *T* = 298 K and *P* = 1 atm. The
temperature and pressure in all simulations were maintained by using
the Nosé–Hoover chain thermostat [60] and the modified Parrinello–
Rahman barostat [61–64]. The quality of equilibration was controlled
both “visually” and statistically: various characteristics of the simulated
models were checked for their stability during the simulation (e.g. unit
cell parameters, density, kinetic and potential components of the total
energy, temperature). After equilibration, a classical MD simulation run
was performed for 1 ns in the *NVE* ensemble to obtain the equilibrium
MD trajectory for further statistical analysis.

Similarly, all mineral-water interfacial models were equilibrated by
classical MD simulations in the *NPT* ensemble for 2 ns at ambient con-
ditions. Then, a classical equilibrium MD simulation run was performed
for 2 ns in the *NVT* ensemble to collect atomic trajectories and calculate
the energies of the systems.

2.3. PIMD simulations details

For the ClayFF-MOH/harmonic model we have also estimated the
nuclear quantum effects using the Centroid Molecular Dynamics (CMD)
approach, which is one of the path integral computational realiza-
tions [65]. This technique enables the extraction of quantum corre-
lation functions using solely the centroid bead. In this method, the
first (centroid) bead, which is a center of the quantum “polymer”, is
uncoupled from the thermostat [65].

The CMD simulations for the system described above were per-
formed with the LAMMPS package using the fix CMD command. The
number of Nosé–Hoover thermostat chains [66] was 6 for all runs. We
considered 16 beads in all of the CMD simulations. The integration
timestep was 0.01 fs and the total number of simulation steps was
100,000, which corresponds to 100 ps. The CMD calculations were
performed only for the harmonic version of ClayFF-MOH, since the
Morse version was shown to be numerically unstable.

A more detailed description of the PIMD method is provided in the
Electronic Supplementary Information.

2.4. DFT simulations details

DFT calculations were performed using the Gaussian and plane
wave basis approach, as implemented in the CP2K simulation software
package [67]. The core electrons were described using the Göedecker–
Teter–Hütter (GTH) [68] pseudopotentials combined with the double-
zeta valence polarized (DZVP MOLOPT) basis set and the auxiliary
plane wave basis set with cutoff energy of 600 Ry. The generalized
gradient approximation (GGA) parameterized by Perdew, Burke, and
Ernzerhof (PBE) [69] was used for the exchange–correlation terms with
the Grimme D3 dispersion correction [70].

First, the geometry and cell optimization was performed using
the Broyden–Fletcher–Goldfarb–Shanno method [71]. After this step,
the model was equilibrated with QMD simulations within the Born–
Oppenheimer approximation in the *NPT* ensemble for 1 ps at *T* = 298
K and *P* = 1 atm. The equations of motion were integrated with the 0.5
fs timestep. Following the equilibration, the *NVE* simulation was run for
25.0 ps, and the atomic coordinates were recorded every timestep for
further statistical analysis.

**3. Results and discussion**

3.1. Crystallographic parameters

The calculated crystallographic parameters of thaumasite from the
classical and quantum MD simulations are shown in Table 1 together
with available experimental data and previous modeling results [1,2,
5,23,25,26]. To obtain these results, classical MD simulations with the
ClayFF/Morse, ClayFF-MOH/harmonic, and ClayFF-MOH/Morse were
performed in the *NPT* ensemble without any symmetry constraints, and
all unit cell parameters were allowed to vary. ClayFF/Morse model has
shown to be unstable, as can be seen in Fig. S4 from the Supplementary
Material. However, both the harmonic and Morse versions of ClayFF-
MOH provide stable structures over the entire simulation time. This
confirms the importance of the MOH term for the stability of hydroxy-
lated crystal structures. Both MOH models yield similar results that are
close to the experimental observations. The discrepancies of the unit
cell parameters are less than 5%. Moreover, our ClayFF-MOH results
show better agreement with the experimental data compared to the
results of the ReaxFF calculations [25].

Our QMD results also agree well with the experimental data. Al-
though all cell parameters were allowed to vary in the *NPT*-QMD
simulations, the deviations of the cell angles are very minor. The recent
DFT calculations [1,26] also agree well with the experimental data,
while the earliest research by Drábik et al. [23] shows substantial
discrepancies, probably due to the different choice of pseudopotentials
and/or exchange-correlation functional.

3.2. Structural parameters

The calculated radial distribution functions (RDF) for the Si–O<sub>h</sub> and Ca–O (both O<sub>w</sub> and O<sub>h</sub>) pairs are shown in Fig. 1. According to experimental data, the average distance between the Si and O<sub>h</sub> atoms is 1.78 Å [2,5]. The calculated average Si–O<sub>h</sub> distance is 1.81 and 1.90 Å for quantum and classical MD runs, respectively. The deviations from the experiment are less than 7%. However, the classical MD model is able to consistently reproduce silicon in six-fold coordination with the hydroxyl groups (Fig. S5).

Calcium has eight oxygen atoms in the first coordination shell of
the thaumasite structure, four of which are O<sub>w</sub> (H<sub>2</sub>O molecules), and
four O<sub>h</sub> (hydroxyls) [5]. The average experimental distance between Ca
and O (Ca–O<sub>w</sub> and Ca–O<sub>h</sub>) is 2.47 Å [2,5]. Our QMD simulations show
good agreement with the experiment, while the classical MD slightly
overestimates the observed distance (2.55 Å), as shown in Fig. 1(B).

According to Drábik et al. [23], thaumasite has four donor–acceptor
pairs of H-bonds: structural hydroxyls – H<sub>2</sub>O molecules (H<sub>o</sub>···O<sub>w</sub> pair);
H<sub>2</sub>O molecules – CO<sub>3</sub><sup>2-</sup> and SO<sub>4</sub><sup>2-</sup> anions (H<sub>w</sub>···O<sub>c</sub> and H<sub>w</sub>···O<sub>s</sub> pairs,
respectively); H<sub>2</sub>O – H<sub>2</sub>O molecules (H<sub>w</sub>···O<sub>w</sub> pair).

Fig. 2 shows the calculated RDFs for H<sub>o</sub>–O<sub>w</sub>, H<sub>w</sub>–O<sub>c</sub>, and H<sub>w</sub>–O<sub>s</sub>
pairs. In all cases, the first peaks from the QMD results are located at
larger distances than the classical MD results. Overall, the discrepancy
in the average H–O distances between quantum and classical MD results
is less than 10%. The PIMD method slightly broadens the peak widths
and shifts them to higher distances. This is especially noticeable for the
H<sub>w</sub>–O<sub>c</sub> and H<sub>w</sub>–O<sub>s</sub> pairs, as shown in Fig. 2(B) and (C).

Similarly to Drábik et al. [23], we can assume that the H<sub>w</sub>···O<sub>c</sub>
and H<sub>w</sub>···O<sub>s</sub> H-bonds are stronger than the H<sub>o</sub>···O<sub>w</sub> bonds because of
their shorter distances. We also find that the H-bonds between H<sub>2</sub>O
molecules (H<sub>w</sub>···O<sub>w</sub> pairs) are quite weak and do not play a crucial
role in the structure of thaumasite.
<table><caption>**Table 1**
Crystallographic unit cell parameters and density of thaumasite: 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><th>Density, g/cm<sup>3</sup></th></tr></thead><tbody><tr><td>ClayFF-MOH/harmonic<sup>a</sup></td><td>10.797</td><td>10.791</td><td>10.879</td><td>1097.40</td><td>89.99</td><td>90.04</td><td>120.02</td><td>1.884</td></tr><tr><td>ClayFF-MOH/Morse<sup>a</sup></td><td>10.776</td><td>10.784</td><td>10.928</td><td>1099.74</td><td>90.03</td><td>89.97</td><td>120.00</td><td>1.880</td></tr><tr><td>DFT (PBE+D3)<sup>a</sup></td><td>10.995</td><td>10.991</td><td>10.457</td><td>1094.46</td><td>90.00</td><td>90.00</td><td>120.00</td><td>1.889</td></tr><tr><td>QMD (PBE+D3)<sup>a</sup></td><td>11.092</td><td>11.078</td><td>10.521</td><td>1116.67</td><td>89.87</td><td>90.14</td><td>120.68</td><td>1.852</td></tr><tr><td>ReaxFF [25]</td><td>11.4048</td><td>11.4048</td><td>10.8245</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>DFT (USP)<sup>b</sup> [23]</td><td>10.89</td><td>10.89</td><td>10.25</td><td>1053</td><td>90</td><td>90</td><td>120</td><td>1.964<sup>a</sup></td></tr><tr><td>DFT (PBE) [1]</td><td>11.071</td><td>11.071</td><td>10.506</td><td>1115.17</td><td>90</td><td>90</td><td>120</td><td>1.854<sup>a</sup></td></tr><tr><td>DFT (PBE) [26]</td><td>11.06</td><td>11.06</td><td>10.52</td><td>1115.11</td><td>90</td><td>90</td><td>120</td><td>1.854<sup>a</sup></td></tr><tr><td>X-ray [5]</td><td>11.0538</td><td>11.0538</td><td>10.4111</td><td>1101.67</td><td>90</td><td>90</td><td>120</td><td>1.876</td></tr><tr><td>Neutron diffraction [2]</td><td>11.0545</td><td>11.0545</td><td>10.4131</td><td>1102.02<sup>a</sup></td><td>90</td><td>90</td><td>120</td><td>1.88</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Calculated radial distribution functions: Si–O<sub>n</sub> pair (A), Ca–O (both O<sub>w</sub> and O<sub>n</sub>) pair (B). Vertical dashed line represents the experimentally determined average interatomic distances [2,5]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 3.3. Dynamics of H-bonding

The H-bond populations for various donor–acceptors pairs (H<sub>w</sub>···O<sub>c</sub>, H<sub>w</sub>···O<sub>s</sub>, H<sub>o</sub>···O<sub>w</sub>) were determined according to the common geometric definition [72]. The continuous-time correlation functions of the H-bonds, C<sub>HB</sub>, were calculated to estimate their lifetimes τ. The same method has already been applied to ettringite [32]. The resulting C<sub>HB</sub> functions are shown in Fig. 3. Table 2 gives the values of the H-bond lifetimes for classical MD and QMD runs, which were obtained by integrating the correlation functions.

Different computational approaches show inconsistencies of the H-bond lifetimes for some donor–acceptor pairs. For example, H<sub>w</sub>···O<sub>s</sub> bonds seems to be stronger than H<sub>w</sub>···O<sub>c</sub> for the QMD model and vice versa for the classical models. Nevertheless, both classical and quantum MD results show that the H-bonds between H<sub>2</sub>O molecules and

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Calculated radial distribution functions for pairs of donors and acceptors of H-bonds: H<sub>o</sub>–O<sub>a</sub> pair (A), H<sub>w</sub>–O<sub>c</sub> pair (B), H<sub>w</sub>–O<sub>s</sub> pair (C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<footnote>a Calculated in this work.</footnote>

<footnote>b Ultrasoft pseudopotentials (USP).</footnote>
<table><caption>**Table 2**
Calculated average lifetimes of intercrystalline H-bonds in thaumasite.</caption><thead><tr><th>Donor-acceptor pair</th><th>Model</th><th>Time (ps)</th></tr></thead><tbody><tr><td rowspan="4">H<sub>w</sub>···O<sub>c</sub></td><td>ClayFF-MOH/Morse</td><td>2.89</td></tr><tr><td>ClayFF-MOH/harmonic</td><td>3.29</td></tr><tr><td>ClayFF-MOH/harmonic+PIMD</td><td>1.95</td></tr><tr><td>QMD (PBE+D3)</td><td>3.37</td></tr><tr><td rowspan="4">H<sub>w</sub>···O<sub>s</sub></td><td>ClayFF-MOH/Morse</td><td>5.21</td></tr><tr><td>ClayFF-MOH/harmonic</td><td>5.22</td></tr><tr><td>ClayFF-MOH/harmonic+PIMD</td><td>3.82</td></tr><tr><td>QMD (PBE+D3)</td><td>2.28</td></tr><tr><td rowspan="4">H<sub>v</sub>···O<sub>w</sub></td><td>ClayFF-MOH/Morse</td><td>0.34</td></tr><tr><td>ClayFF-MOH/harmonic</td><td>0.55</td></tr><tr><td>ClayFF-MOH/harmonic+PIMD</td><td>0.29</td></tr><tr><td>QMD (PBE+D3)</td><td>0.27</td></tr></tbody></table>

anions (H<sub>w</sub>···O<sub>c</sub> and H<sub>w</sub>···O<sub>s</sub> pairs) have longer lifetimes than the H-bonds between hydroxyls and the H<sub>2</sub>O molecules (H<sub>v</sub>···O<sub>w</sub> pair). This confirms the argument above that the H-bonds between H<sub>2</sub>O molecules and anions are stronger. The modeling of ettringite [32] has also shown that the H-bonds between structural hydroxyls and H<sub>2</sub>O molecules are weaker than the H-bonds between the H<sub>2</sub>O molecules and the SO<sub>4</sub><sup>2-</sup> anions. However, ettringite exhibits a denser network of H-bonds due to the greater number of water molecules in its structure [73]. Thus, it is likely that the stability of ettringite when exposed to aqueous salt solutions should be higher than that of thaumasite in sulfate attack scenarios.

A comparison of the PIMD and classical MD results shows that the correlation functions decay faster (Fig. 3) and that the lifetimes of H-bonds are shorter for certain donor-acceptor pairs. This is also observed for clay minerals (kaolinite and pyrophyllite) [44]. This indicates a weakening of the H-bonding network as a result of hydrogen delocalization. In general, our present results are consistent with the recent observation that taking into account the nuclear quantum effects in water is significant for the H-bond dynamics but not as important for the average structure [74].

### 3.4. Vibrational properties

The so-called power spectra (PS) of different modes of classical atomic vibrations in the crystal provide another useful way to assess the contribution of H-bonding to the properties of thaumasite. The PS were obtained by the Fourier transform of the velocity autocorrelation function of the respective atoms [32].

Fig. 4(A) shows the total PS of thaumasite from quantum and classical MD simulations. The peaks reflecting the O–H stretching modes of hydroxyls and H<sub>2</sub>O molecules are located in the high-frequency region of the total PS. Therefore, decomposed PS were calculated to better understand different vibrational contributions of H<sub>2</sub>O molecules and hydroxyls to the total spectrum (Fig. 4(B) and (C), respectively).

Experimental Fourier Transform Infrared (FTIR) spectra [1,2] show two peaks that reflect the symmetric and asymmetric stretching vibrations of the O–H bonds for both H<sub>2</sub>O molecules and hydroxyls. The modes at 3370 and 3460 cm<sup>-1</sup> correspond to symmetric and asymmetric stretching of the H<sub>2</sub>O molecules (Fig. 4(B)), respectively [2]. The best agreement with the experimental data is achieved using the QMD simulation, which gives the peak maximum at 3431 cm<sup>-1</sup>. All peaks obtained in classical MD simulations are red-shifted, which is a known problem with the SPC and SPC/E water models. PIMD methodology slightly improves the agreement with the FTIR data, resulting in the peak maximum of this mode at 3640 cm<sup>-1</sup>.

The spectral peaks obtained in PIMD simulations are weaker and broader compared to those obtained in classical MD simulations with the ClayFF-MOH/harmonic model, which is the result of hydrogen delocalization. The Morse potential also affords peak positions slightly

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

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

<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Calculated continuous time autocorrelation functions of intracrystalline H-bonds in thaumasite: H<sub>w</sub>···O<sub>c</sub> donor-acceptor pair (A), H<sub>w</sub>···O<sub>s</sub> donor-acceptor pair (B), H<sub>v</sub>···O<sub>w</sub> donor-acceptor pair (C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

closer to the experimental values. However, classical models do not achieve the same level of accuracy as the QMD results.

The symmetric and asymmetric stretching vibrational modes of hydroxyls (Fig. 4(C)) are located at 3340 and 3500 cm<sup>-1</sup>, respectively [2]. These peaks occur at 3373 and 3473 cm<sup>-1</sup> in the QMD and ClayFF-MOH/Morse models, showing a good agreement with experimental data. The ClayFF-MOH/harmonic model yields the most shifted peak, which again can be slightly improved with the addition of the PIMD model. Again, the PIMD peak of hydroxyls are weaker and broader. Similar behavior of the O–H peaks was also noted for clays [43,44].

For all types of simulations, the peaks reflecting the H<sub>2</sub>O angular bending (deformation) mode at ~1651 cm<sup>-1</sup> (Fig. 4(A)) are in excellent agreement with the FTIR data [1,2]. The peaks reflecting the O–C and
<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Power spectra of thaumasite: total spectrum (A), vibrations of H₂O molecules in the high-frequency region (B), vibrations of hydroxyls (OH) in the high-frequency region (C). Vertical dashed lines indicate the experimentally determined symmetric and asymmetric O–H stretching frequencies [2]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

O–S stretching modes of the CO₂⁻ and SO₄²⁻ anions are located at ~1397 and ~1067 cm⁻¹, respectively (Fig. S6).

### 3.5. Mechanical properties of thaumasite

The full elastic stiffness tensor of the thaumasite crystal was calculated from the series of separate equilibrium classical MD trajectories with the ClayFF-MOH/Morse and ClayFF-MOH/harmonic models, similarly to the earlier ettringite simulations [32]. The bulk modulus (K), shear modulus (G), Young's modulus (E), and the Poisson ratio (μ) of thaumasite were also obtained using the Voigt–Reuss–Hill (VRH)

approximation [75]. The comparison of our results with available experimental data and previous ReaxFF and DFT simulations is presented in Table 3.

The calculated mechanical properties with the ClayFF-MOH/Morse and ClayFF-MOH/harmonic models are close to each other. However, the discrepancy between the calculated C₁₁, C₁₂, and C₃₃ elastic constants and the Brillouin spectroscopic data [24] is up to 32%. Nevertheless, other calculated constants agree well with the experiment. The underestimation of the elastic constants is reflected in the calculated G and E moduli, which also show deviations of up to 25%. However, the calculated bulk modulus and Poisson ratio are close to the experimental values.

Previous classical MD simulations using the ReaxFF force field [25] have also underestimated some mechanical characteristics of thaumasite. The maximum discrepancy with the experimental data was up to 58% for the C₆₆ constant [25]. In contrast, the DFT calculations of Scholtzová et al. [24] and San et al. [26] overestimate the mechanical characteristics of thaumasite. The maximum overestimations reach from 30% to 44% for the C₆₆ constant obtained from the former and latter calculations, respectively. However, the average elastic constant overestimations from these calculations are between 12% and 20%, respectively. Our calculations with the empirical ClayFF-MOH model result in the average elastic constant underestimations of 23% and 17% for simulations with Morse and harmonic potential, respectively. Thus, DFT calculations are only slightly better in reproducing the mechanical characteristics of thaumasite, and one can conclude that the ClayFF-MOH/harmonic model offers a rather satisfactory representation of the thaumasite elastic properties.

### 3.6. Thermodynamic properties

Helmholtz free energy (F), specific heat at constant volume (C_v), and entropy (S) of thaumasite (Fig. 5) were obtained from the phonon frequencies using the Phonopy code [76]. The initial crystal structure was prepared by geometry optimization under the proper symmetry constraints. The uniform displacement of 0.01 Å was set for the finite difference method. For the classical simulation, the 4 × 4 × 4 supercell was selected, while for the DFT calculations only the primitive cell was used. The small cell size for DFT calculations is an acceptable choice, as it has already been shown for ettringite [77]. Fig. 5 demonstrates that the results for both the ClayFF-MOH/Morse and DFT simulations of thaumasite are very close to each other for all thermodynamic quantities calculated over a wide range of temperatures. Previous calculations of the thermodynamic quantities for ettringite have shown that the results of classical modeling with the original ClayFF parametrization (without the MOH term) agree well with the DFT results [77]. Assuming that the presently used ClayFF modification slightly affects the thermodynamics of thaumasite, we present here only the results for the ClayFF-MOH/Morse model.

To compare the Helmholtz free energy for both crystals, we also set the Helmholtz free energy value at 0 K for thaumasite as for ettringite [77]. As a result, ettringite is seen as more stable than thaumasite (Fig. 5(A)). This is consistent with the thermodynamic database for hydrated Portland cements [78,79], where the Gibbs energy of ettringite (Δ<sub>f</sub> G° = -15205.94 kJ/mol) is lower than that of thaumasite (Δ<sub>f</sub> G° = -7564.52 kJ/mol).

Not surprisingly, the calculated entropy of thaumasite at room temperature does not agree well with the experimental value (see Table S1). This is also observed for ettringite (Fig. 5(C)). The lower entropy of thaumasite indicates that it is more ordered than ettringite, which is also associated with the higher crystal symmetry of thaumasite. A comparison of the calculated heat capacity with experimental value also shows a deviation for thaumasite (see Table S1). However, for both crystals, the discrepancies in the calculated entropy and heat capacity with the experimental values are less than 43% [3,78,79].
<table><caption>**Table 3**

Mechanical properties of thaumasite: simulation results and experimental data.</caption><thead><tr><th></th><th>ClayFF-MOH/Morse<sup>b</sup></th><th>ClayFF-MOH/harmonic<sup>b</sup></th><th>ReaxFF [25]</th><th>DFT(PBE) [26]</th><th>DFT(PBE) [24]</th><th>Experiment [24]</th></tr></thead><tbody><tr><td>C<sub>11</sub><sup>a</sup></td><td>47.7</td><td>47.0</td><td>41.26</td><td>68.10</td><td>67.3</td><td>59.9 ± 0.1</td></tr><tr><td>C<sub>12</sub><sup>a</sup></td><td>25.8</td><td>28.1</td><td>32.88</td><td>31.83</td><td>34.1</td><td>34.3 ± 0.2</td></tr><tr><td>C<sub>13</sub><sup>a</sup></td><td>21.0</td><td>20.4</td><td>25.98</td><td>24.70</td><td>24.2</td><td>24.1 ± 0.3</td></tr><tr><td>C<sub>33</sub><sup>a</sup></td><td>41.9</td><td>42.9</td><td>65.84</td><td>67.07</td><td>66.6</td><td>61.4 ± 0.3</td></tr><tr><td>C<sub>44</sub><sup>a</sup></td><td>10.4</td><td>12.3</td><td>8.77</td><td>19.57</td><td>16.8</td><td>13.9 ± 0.2</td></tr><tr><td>C<sub>66</sub><sup>a</sup></td><td>10.2</td><td>11.9</td><td>5.44</td><td>18.44</td><td>16.6</td><td>12.8 ± 0.2</td></tr><tr><td>K<sup>a</sup></td><td>36.9</td><td>36.6</td><td>35.33</td><td>40.45</td><td>40.6</td><td>38.5 ± 0.3</td></tr><tr><td>G<sup>a</sup></td><td>10.9</td><td>10.9</td><td>8.31</td><td>19.53</td><td>17.8</td><td>14.5 ± 0.4</td></tr><tr><td>E<sup>a</sup></td><td>29.7</td><td>29.3</td><td>23.12</td><td>50.46</td><td>46.1</td><td>38.6 ± 0.8</td></tr><tr><td>ν</td><td>0.36</td><td>0.36</td><td>0.39</td><td>0.29</td><td>0.31</td><td>0.33 ± 0.01</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Thaumasite thermodynamic properties: (A) Helmholtz free energy, (B) Heat capacity at constant volume, (C) Entropy. Dashed green line - DFT and ClayFF simulations for ettringite [77]. Green square and diamond - experimental $C_p$ and S for ettringite [79], respectively. Blue square and diamond - experimental $C_p$ [78] and S [3] for thaumasite, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 3.7. Thaumasite and ettringite interfaces with water

### 3.7.1. Adsorption of water molecules

The structure of water at different crystal faces of thaumasite and ettringite was investigated by classical MD simulations using the ClayFF-MOH/harmonic model. For the Morse version, a smaller timestep was required to avoid numerical instabilities. Taking into account that both versions of the ClayFF-MOH force field show very similar results for the crystal properties of thaumasite, we expect that they should work equally well for the simulations of mineral-water interfaces. Therefore, only the harmonic version was used in these calculations.

The interaction of aqueous salt solutions with ettringite surfaces has already been studied before [32,80–82]. Ettringite needles exhibit six large elongate (1 0 0) and (0 1 0) rectangular prism faces and two flat hexagonal (0 0 1) end faces [83]. The (0 1 0) and (1 0 0) surfaces are identical due to the crystal space symmetry of ettringite [84]. The structural properties of the interfacial aqueous solution, the diffusional mobility of the surface ions and H₂O molecules, and the corresponding adsorption isotherms were mainly studied for the ettringite (1 0 0) surface by classical MD simulations. However, the (0 0 1) surface is also important and has been recently investigated using the DFT methodology [84].

In our simulations, we initially focus on the (1 0 0) interface. The profiles of atomic densities for $H_ω$ and $O_ω$ of H₂O molecules as functions of distance from the crystal surface are shown in Fig. 6. The plane defined by the time-averaged positions of the sulfate oxygen atoms ($O_s$) was taken as a reference boundary of the mineral surface (vertical dashed line (see also Fig. S7)). Water adsorption occurs on the surfaces of both ettringite and thaumasite, with noticeable oscillations in water density near these surfaces. This is indicated by the intense $O_ω$ peaks at 1.9 and 1.5 Å for ettringite and thaumasite, respectively. H₂O molecules from solution are adsorbed mainly on the columns of thaumasite and ettringite due to the acceptance of H-bonds from the surfaces. H₂O molecules are also adsorbed on the sulfate and carbonate anions by donating H-bonds to the oxygen atoms of the anions. The $H_ω$ peaks located at 1.0 and 0.8 Å for ettringite and thaumasite are indicative of this.

The aforementioned peaks of the $H_ω$ and $O_ω$ interfacial profiles for thaumasite are more intense and narrower than those for ettringite. Hence, we can assume that H₂O molecules adsorb more strongly on the thaumasite surface. However, surprisingly, the integral of the water density profile gives the same amount of water for both thaumasite and ettringite.

Fig. 7 shows the one-dimensional density profiles of $H_ω$ and $O_ω$ near the (0 0 1) surfaces of thaumasite and ettringite. The vertical dashed line indicates the reference boundary of the mineral surface, which is determined by the time-averaged positions of the uttermost hydroxyl oxygen atoms $O_ħ$ (see also Fig. S8). As for the (1 0 0) surfaces, the water density profiles are shifted closer to the surface of thaumasite, and the density peaks are more intense and narrower. This indicates a greater localization of H₂O molecules at the surface of thaumasite in

<footnote>a In GPa.</footnote>

<footnote>b Calculated in this work.</footnote>
<table><caption>**Table 4**
Calculated solid-water interfacial energies (Δ*E*<sub>SL</sub>/A) for (h k l) surfaces of thaumasite and ettringite.</caption><thead><tr><th>Mineral</th><th>(h k l)</th><th>ΔE<sub>SL</sub>/A (mJ/m<sup>2</sup>)</th></tr></thead><tbody><tr><td rowspan="2">Thaumasite</td><td>(1 0 0)</td><td>-7 ± 13</td></tr><tr><td>(0 0 1)</td><td>-40 ± 20</td></tr><tr><td rowspan="2">Ettringite</td><td>(1 0 0)</td><td>9 ± 20</td></tr><tr><td>(0 0 1)</td><td>34 ± 23</td></tr></tbody></table>

comparison with ettringite. However, again, the integral of the water density profiles from 0 to 10 Å gives the same number of interfacial H₂O molecules for both thaumasite and ettringite.

Fig. S7 and S8 also show that the H₂O molecules from the interfacial aqueous phase could effectively penetrate into the crystal structure, while the structural H₂O molecules could occasionally diffuse in the opposite direction. This exchange process indicates the dynamical equilibrium of water molecules near the surface.

### 3.7.2. Interfacial free energy

The solid-water interfacial energy Δ*E*<sub>SL</sub>/A was obtained using the three-box method [85,86]. The calculated values for thaumasite and ettringite are given in Table 4. We neglect the entropy of forming the solid-liquid interface Δ*S*<sub>SL</sub> from bulk solid and bulk water. However, for calcium sulfate minerals (anhydrite, bassanite, gypsum), the average loss of entropy was estimated at 15 ± 5 mJ/m² [86].

Both thaumasite (1 0 0) and (0 0 1) surfaces exhibit negative interfacial energy of -7 and -40 mJ/m², respectively. However, both ettringite (1 0 0) and (0 0 1) surfaces have positive solid-water interfacial energy of 9 and 34 mJ/m², respectively. Thus, the affinity of ettringite for water is lower. The entropy contribution Δ*S*<sub>SL</sub> for thaumasite might be slightly larger due to the greater localization of H₂O molecules near surfaces, as shown in Figs. 6 and 7. Nevertheless, the solid-water interfacial tensions for thaumasite are lower than those for ettringite.

The more negative solid-water interfacial energy for the thaumasite (0 0 1) surface indicates that this surface is more soluble. This is probably due to the smaller number of H₂O molecules in the thaumasite structure. Mishra et al. [86] showed that negative solid-liquid interfacial energies also correlate with the presence of voids and water channels in the bassanite structure. Furthermore, the stability of ettringite decreases as the number of structural H₂O molecules decreases [87,88]. Thus, thaumasite is more soluble than ettringite. It is also confirmed by the thermodynamic data database for hydrated Portland cements [79], where the solubility constant of ettringite (log*K*<sub>S0</sub> = -44.90) is lower than that of thaumasite (log*K*<sub>S0</sub> = -24.75).

### 3.8. Thaumasite and ettringite interfaces

At both temperatures 298 and 278 K the thickness of the water adsorption layer on the ettringite (1 0 0) surface is approximately 10 Å (Fig. 6 and Fig. S9, respectively). The diameter of the thaumasite column is also about 9–10 Å. Therefore, the thickness of the adsorption layer of water is sufficient for the nucleation of the thaumasite and formation of its column. Thus, we suggest that epitaxial growth of thaumasite is possible on the ettringite surface.

In epitaxial crystal growth, nucleation is not trivial and depends on many factors (e.g., solution supersaturation, diffusion and mobility of adatoms and admolecules on the surface, etc.). Thus, epitaxial crystal growth requires additional sophisticated series of simulations. However, even based on the results presented above, we can estimate the stability of the mixture of thaumasite and ettringite crystals from classical MD simulations.

An atomistic model of a mixture of thaumasite and ettringite crystals is shown in Fig. 8. The initial lateral dimensions of the supercell were adopted as for ettringite, since its dimensions are slightly larger

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

<figure><img src="image_12.png" /><figcaption>**Fig. 6.** Atomic density profiles of O<sub>w</sub>(A) and H<sub>w</sub>(B) in aqueous solution at the (1 0 0) plane surface of thaumasite and ettringite from classical MD simulations with the ClayFF-MOH force field. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

than those of thaumasite. All simulations details were similar to the ones used for the mineral-water models. Over the entire simulation time, the mixed model remained stable showing only a very slight restructuring of cell sizes in comparison to the pure crystal (see Table S2). No significant changes in the column structure of thaumasite and ettringite are observed at the interface. The calculated solid-solid interfacial energy at ambient conditions is 46 ± 16 mJ/m². The positive interfacial energy and the absence of significant changes in the structure of both crystals at the interface indicate the possible stable existence of such a crystal mixture [89].

## 4. Conclusions

In this work, we used classical molecular dynamics simulations with the ClayFF-MOH model [29] to study the properties of thaumasite crystals and their aqueous interfaces. The inclusion of explicit Ca–O–H and Si–O–H bending terms in the simulation model stabilizes the crystal structure, which was not observed for the original ClayFF version [27].

The computed structural, mechanical, and thermodynamic properties of thaumasite show good agreement with available experimental and with the results of DFT calculations. As expected, the application of Morse potential for the stretching vibrations of the structural hydroxyls and H₂O molecules allows to improve the reproduction of vibrational properties. Path Integral Molecular Dynamics (PIMD) with the ClayFF-MOH/harmonic version of the force field also improves the accuracy of the calculated power spectra by shifting the O–H stretching modes closer to the experimental frequencies.

Investigation of the intracrystalline hydrogen bonding dynamics shows that the most stable and strongest H-bonds are observed between
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 7. Atomic density profiles of O<sub>w</sub>(A) and H<sub>w</sub>(B) in aqueous solution at the (0 0 1) plane surface of thaumasite and ettringite from classical MD with the ClayFF-MOH force field. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 8. Equilibrium snapshot of an atomistic model of a mixture of thaumasite and ettringite crystals. Atom labels: O (red), H (pale gray), Si (orange), Ca (green), S (yellow), C (dark gray). OVITO software [90] was used for visualization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the structural H₂O molecules and the anions (CO₃²⁻ and SO₄²⁻). The H-bonds between hydroxyls and H₂O molecules in the crystal structure are weaker, as is was previously emphasized by Drábik et al. [23]. The same result was also observed for ettringite [32]. This is consistent with the thermal analysis of the decomposition of thaumasite and ettringite [3,23,91,92]. Breaking of the H-bonds between structural water molecules and Si(OH)₂⁻ or Al(OH)₃⁻ (thaumasite and ettringite respectively) occurs even at relatively low temperatures (363–393 K). This promotes the release of water and collapse of the Si(OH)₂⁻ and Al(OH)₃⁻ octahedra and a rearrangement of the crystal structure of thaumasite and ettringite at high temperatures [92,93].

Additionally, the more realistic calculations with PIMD showed a weakening of the H-bonding network within the thaumasite structure, suggesting that further application of PIMD would be even more interesting in modeling the mineral-aqueous solution interfaces. Obviously, the weakening of the H-bonding network would affect the structural and transport properties of aqueous solutions near mineral surfaces. This, in turn, would affect the stability and solubility of thaumasite and ettringite.

Water structuring at the (1 0 0) and (0 0 1) surfaces for thaumasite and ettringite shows that the affinity of thaumasite for water is greater. However, the total amount of H₂O molecules on the surface is the same for both minerals. The MD simulations also demonstrate that the solid (1 0 0) interface between the thaumasite and ettringite crystals is stable confirming the assumption that epitaxial growth of thaumasite on the ettringite surface is possible [4]. Of course, this observation does not exclude that other ways of thaumasite formation (direct precipitation from solution or formation via solid-solution with ettringite) are also possible. Which pathway of thaumasite formation is more predominant remains a subject of further detailed research.

Present calculations of the (h k l) solid-water interfacial energies confirm that thaumasite is more soluble than ettringite. The (0 0 1) surface of thaumasite is more susceptible to dissolution than its (1 0 0) surface. Due to the smaller number of H₂O molecules in the structure of thaumasite, the density of the H-bonding network is lower at the thaumasite (0 0 1) surface. Dissolution of structural H₂O molecules from the thaumasite (0 0 1) surface is most likely, as intracrystalline H-bonds between structural hydroxyls and H₂O molecules are the weakest. This suggestion is also relevant for ettringite.

### CRediT authorship contribution statement

E.V. Tararushkin: Writing – review & editing, Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Data curation. G.S. Smirnov: Writing – review & editing, Project administration, Formal analysis. N.D. Kondratyuk: Writing – original draft, Investigation. A.G. Kalinichev: Writing – review & editing, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The article was prepared within the framework of the HSE University Basic Research Program. This research was supported through computational resources of HPC facilities at HSE University [94]. N.D.K. is grateful for the support from the Project No. FSMG-2025-0006 (PIMD simulations for PS). 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. We are grateful to Denis Potapov for his assistance with the PS calculations for PIMD.
# Appendix A. Supplementary data

Supplementary material related to this article (simulated model snapshots and additional figures) can be found online at https://doi.org/10.1016/j.cemconres.2025.107944.

## Data availability

Data will be made available on request.

## References

[1] E. Scholtzová, L. Kucková, J. Kožíšek, H. Pálková, M. Tunega, Experimental and computational study of thaumasite structure, Cem. Concr. Res. 59 (2014) 66–72, http://dx.doi.org/10.1016/j.cemconres.2014.02.002.

[2] G.D. Gatta, G.J. McIntyre, J.G. Swanson, S.D. Jacobsen, Minerals in cement chemistry: A single-crystal neutron diffraction and Raman spectroscopic study of thaumasite, Ca₃Si(OH)₆(CO₃)(SO₄)·12H₂O, Am. Mineral. 97 (7) (2012) 1060–1069, http://dx.doi.org/10.2138/am.2012.4058.

[3] Y.D. Gritsenko, M.F. Vigasina, L.V. Mel'chakova, L.P. Ogorodova, D.A. Ksenofontova, S.K. Dedushenko, Thermal and thermochemical study of thaumasite, Geochem. Int. 61 (12) (2023) 1273–1282, http://dx.doi.org/10.1134/S0016702923110046.

[4] S. Köhler, D. Heinz, L. Urbanas, Effect of ettringite on thaumasite formation, Cem. Concr. Res. 36 (2006) 697–706, http://dx.doi.org/10.1016/j.cemconres.2005.11.006.

[5] S.D. Jacobsen, J.R. Smyth, R.J. Swope, Thermal expansion of hydrated six-coordinated silicon in thaumasite, Ca₃Si(OH)₆(CO₃)(SO₄)·12H₂O, Phys. Chem. Miner. 30 (2003) 321–329, http://dx.doi.org/10.1007/s00269-003-0328-0.

[6] R.A. Edge, H.F.W. Taylor, Crystal structure of thaumasite, [Ca₃Si(OH)₆·12H₂O](SO₄)(CO₃), Acta Crystallogr. Sect. B: Struct. Sci. Cryst. Eng. Mater. 27 (1971) 594–601, http://dx.doi.org/10.1107/S0567740871002619.

[7] J.A. Duffy, D.E. Macphee, The coordination number of silicon in silicon-oxygen compounds: The special case of 6-fold coordination in thaumasite, J. Phys. Chem. B 111 (30) (2007) 8740–8745, http://dx.doi.org/10.1021/jp071343n.

[8] Q. Williams, R.J. Hemley, M.B. Kruger, R. Jeanloz, High-pressure infrared spectra of α-quartz, coesite, stishovite and silica glass, J. Geophys. Res. Solid Earth 98 (B12) (1993) 22157–22170, http://dx.doi.org/10.1029/93JB02171.

[9] G. Gibbs, A. Wallace, D. Cox, R. Downs, N. Ross, K. Rosso, Bonded interactions in silica polymorphs, silicates, and siloxane molecules, Am. Mineral. 94 (8–9) (2009) 1085–1102, http://dx.doi.org/10.2138/am.2009.3215.

[10] J.F. Stebbins, S. Bista, Pentacoordinated and hexacoordinated silicon cations in a potassium silicate glass: Effects of pressure and temperature, J. Non Cryst. Solids 505 (2019) 234–240, http://dx.doi.org/10.1016/j.jnoncrysol.2018.11.001.

[11] T. Mahadevan, J. Du, Silica, Silicate, and Aluminosilicate Glasses, Wiley, 2022, pp. 186–223, http://dx.doi.org/10.1002/9781118939079.ch7, Ch. 7.

[12] R. Schliemann, S.V. Churakov, Pyrophyllite dissolution at elevated pressure conditions: An ab initio study, Geochim. Cosmochim. Acta 307 (2021) 42–55, http://dx.doi.org/10.1016/j.gca.2021.05.017.

[13] H. Taylor, Cement Chemistry, second ed., Thomas Tedford Publishing, London, UK, 1997.

[14] J. Bensted, Thaumasite–direct, woodfordite and other possible formation routes, Cem. Concr. Compos. 25 (2003) 873–877, http://dx.doi.org/10.1016/S0958-9465(03)00115-X.

[15] F. Bellmann, J. Stark, Prevention of thaumasite formation in concrete exposed to sulphate attack, Cem. Concr. Res. 37 (2007) 1215–1222, http://dx.doi.org/10.1016/j.cemconres.2007.04.007.

[16] F. Bellmann, J. Stark, The role of calcium hydroxide in the formation of thaumasite, Cem. Concr. Res. 38 (2008) 1154–1161, http://dx.doi.org/10.1016/j.cemconres.2008.04.005.

[17] K. Sotiriadis, P. Mácová, A.S. Mazur, A. Viani, P.M. Tolstoy, S. Tsivilis, Long-term thaumasite sulfate attack on portland-limestone cement concrete: A multi-technique analytical approach for assessing phase assemblage, Cem. Concr. Res. 130 (2020) 105995, http://dx.doi.org/10.1016/j.cemconres.2020.105995.

[18] N. Crammond, The occurrence of thaumasite in modern construction – a review, Cem. Concr. Compos. 24 (3–4) (2002) 393–402, http://dx.doi.org/10.1016/S0958-9465(01)00092-0.

[19] J. Aguilera, M. Varela, T. Vázquez, Procedure of synthesis of thaumasite, Cem. Concr. Res. 31 (8) (2001) 1163–1168, http://dx.doi.org/10.1016/S0008-8846(01)00536-1.

[20] T. Schmidt, B. Lothenbach, M. Romer, K. Scrivener, D. Rentsch, R. Figi, A thermodynamic and experimental study of the conditions of thaumasite formation, Cem. Concr. Res. 38 (2008) 337–349, http://dx.doi.org/10.1016/j.cemconres.2007.11.003.

[21] K.L. Scrivener, R.J. Kirkpatrick, Innovation in use and research on cementitious material, Cem. Concr. Res. 38 (2) (2008) 128–136, http://dx.doi.org/10.1016/j.cemconres.2007.09.025.

[22] R. Mishra, A. Mohamed, D. Geissbuhler, H. Manzano, T. Jamil, R. Shahsavari, A. Kalinichev, S. Galmarini, L. Tao, H. Heinz, R. Pellenq, A. van Duin, S. Parker, R. Flatt, P. Bowen, Cemff: A force field database for cementitious materials including validations, applications and opportunities, Cem. Concr. Res. 102 (2017) 68–89, http://dx.doi.org/10.1016/j.cemconres.2017.09.003.

[23] M. Drábik, M. Tunega, S. Balkovic, V.Š. Fajnor, Computer simulations of hydrogen bonds for better understanding of the data of thermal analysis of thaumasite, J. Therm. Anal. Calorim. 85 (2006) 469–475, http://dx.doi.org/10.1007/s10973-005-7236-0.

[24] E. Scholtzová, M. Tunega, S. Speziale, Mechanical properties of ettringite and thaumasite—DFT and experimental study, Cem. Concr. Res. 77 (2015) 9–15, http://dx.doi.org/10.1016/j.cemconres.2015.06.008.

[25] S. Hajilar, B. Shafei, T. Cheng, A. Jaramillo-Botero, Reactive molecular dynamics simulations to understand mechanical response of thaumasite under temperature and strain rate effects, J. Phys. Chem. A 121 (24) (2017) 4688–4697, http://dx.doi.org/10.1021/acs.jpc.a7b02824.

[26] S. San, N. Li, Y. Tao, W. Zhang, W.-Y. C, Understanding the atomic and electronic origin of mechanical property in thaumasite and ettringite mineral crystals, J. Am. Ceram. Soc. 101 (11) (2018) 5177–5187, http://dx.doi.org/10.1111/jace.15774.

[27] R. Cygan, J.-J. Liang, A. Kalinichev, Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field, J. Phys. Chem. B 108 (4) (2004) 1255–1266, http://dx.doi.org/10.1021/jp0363287.

[28] A.G. Kalinichev, J. Wang, R.J. Kirkpatrick, Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: Application to cement materials, Cem. Concr. Res. 37 (2007) 337–347, http://dx.doi.org/10.1016/j.cemconres.2006.07.004.

[29] R. Cygan, J. Greathouse, A. Kalinichev, Advances in Clayff molecular simulation of layered and nanoporous materials and their aqueous interfaces, J. Phys. Chem. C 125 (2021) 17573–17589, http://dx.doi.org/10.1021/acs.jpc.c1c04600.

[30] M. Pouvreau, J. Greathouse, R. Cygan, A. Kalinichev, 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 121 (27) (2017) 14757–14771, http://dx.doi.org/10.1021/acs.jpc.b705362.

[31] M. Pouvreau, J. Greathouse, R. Cygan, A. Kalinichev, 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 123 (18) (2019) 11628–11638, http://dx.doi.org/10.1021/acs.jpcb.9b00514.

[32] E. Tararushkin, V. Pisarev, A. Kalinichev, Atomistic simulations of ettringite and its aqueous interfaces: Structure and properties revisited with the modified ClayFF force field, Cem. Concr. Res. 156 (2022) 106759, http://dx.doi.org/10.1016/j.cemconres.2022.106759.

[33] T. Ho, J. Greathouse, Y. Wang, L. Criscenti, Atomistic structure of mineral nano-aggregates from simulated compaction and dewatering, Sci. Rep. 7 (2017) 15286, http://dx.doi.org/10.1016/j.cemconres.2020.105995.

[34] E. Tararushkin, V. Pisarev, A. Kalinichev, Interaction of nitrite ions with hydrated portlandite surfaces: Atomistic computer simulation study, Materials 16 (14) (2023) 5026, http://dx.doi.org/10.3390/ma16145026.

[35] D. Marx, J. Hutter, AB Initio Molecular Dynamics: Basic Theory and Advanced Methods, CAMBRIDGE UNIVERSITY PRESS, New York, 2009.

[36] T.E. Markland, M. Ceriotti, Nuclear quantum effects enter the mainstream, Nat. Rev. Chem. 2 (3) (2018) 0109, http://dx.doi.org/10.1038/s41570-017-0109-9.

[37] M.E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, 2023, http://dx.doi.org/10.1093/oso/9780198825562.001.0001.

[38] M. Ceriotti, W. Fang, P.G. Kusalik, R.H. McKenzie, A. Michaelides, M.A. Morales, T.E. Markland, Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges, Chem. Rev. 116 (13) (2016) 7529–7550, http://dx.doi.org/10.1021/acs.chemrev.5b00674.

[39] D.O. Potapov, N.D. Kondratyuk, Influence of nuclear quantum effects on the vibrational spectrum and structure of H₂O@C₆₀ endofullerene, JETP Lett. 119 (12) (2024) 923–928, http://dx.doi.org/10.1134/S0021364024601325.

[40] I.S. Novikov, E.M. Makarov, Y.V. Suleimanov, A.V. Shapeev, Towards reliable calculations of thermal rate constants: Ring polymer molecular dynamics for the OH+HBr→Br+ H₂O reaction, Chem. Phys. Lett. (2024) 141620, http://dx.doi.org/10.1016/j.cplett.2024.141620.

[41] C. Liang, A. Rayabharam, N. Aluru, Structural and dynamical properties of H₂O and D₂O under confinement, J. Phys. Chem. B 127 (29) (2023) 6532–6542, http://dx.doi.org/10.1021/acs.jpcb.3c02868.

[42] Q. Li, X. Lu, M. Chen, L. Zhang, Y. Cheng, X. Liu, Z. Yin, Lithium partition and isotopic fractionation for cation exchange in clay: Insights from molecular simulations, Geochim. Cosmochim. Acta 358 (2023) 148–161, http://dx.doi.org/10.1016/j.gca.2023.07.024.

[43] C. Gao, Y. Zhang, Q. Liu, Y. Yang, Y. Liu, Path-integral molecular dynamics predictions of equilibrium h and o isotope fractionations between brucite and water, Geochim. Cosmochim. Acta 346 (2023) 207–222, http://dx.doi.org/10.1016/j.gca.2023.02.004.

[44] P.K. Kurapothula, S. Shepherd, D.M. Wilkins, Hydrogen-bonding and nuclear quantum effects in clays, J. Chem. Phys. 156 (8) (2022) 084702, http://dx.doi.org/10.1063/5.0083075.
[45] P.K. Kurapothula, S. Shepherd, D.M. Wilkins, Competing nuclear quantum effects and hydrogen-bond jumps in hydrated kaolinite, J. Phys. Chem. Lett. 14 (6) (2023) 1542-1547, http://dx.doi.org/10.1021/acs.jpclett.2c03896.

[46] R. Dupuis, J.S. Dolado, M. Benoit, J. Surga, A. Ayuela, Quantum nuclear dynamics of protons within layered hydroxides at high pressure, Sci. Rep. 7 (1) (2017) 4842, http://dx.doi.org/10.1038/s41598-017-04080-2.

[47] S. Schaack, P. Depondt, S. Huppert, F. Finocchi, Quantum driven proton diffusion in brucite-like minerals under high pressure, Sci. Rep. 10 (2020) http://dx.doi.org/10.1038/s41598-020-64813-8.

[48] S. Barnett, D. Macphee, N. Crammond, Extent of immiscibility in the ettringite-thaumasite system, Cem. Concr. Compos. 25 (8) (2003) 851-855, http://dx.doi.org/10.1016/S0958-9465(03)00116-1.

[49] M. Medala, C. Labbez, I. Pochard, A. Nonat, Ettringite surface chemistry: Interplay of electrostatic and ion specificity, J. Colloid Interface Sci. 354 (2011) 765-770, http://dx.doi.org/10.1016/j.jcis.2010.11.031.

[50] T. Missana, O. Almendros-Ginest, F. Colmenero, A.M. Fernández, M. García-Gutiérrez, Investigation of the surface charge behaviour of ettringite: Influence of pH, calcium, and sulphate ions, Heliyon 10 (16) (2024) e36117, http://dx.doi.org/10.1016/j.heliyon.2024.e36117.

[51] A.G. Kalinichev, Atomistic modeling of clays and related nanoporous materials with ClayFF force field, in: Computational Modeling in Clay Mineralogy, in: AIPEA Educational Series, Bari, Italy, 2021, Pub. No. 3.

[52] S. Kerisit, C. Liu, Molecular simulation of the diffusion of uranyl carbonate species in aqueous solution, Geochim. Cosmochim. Acta 74 (17) (2010) 4937-4952, http://dx.doi.org/10.1016/j.gca.2010.06.007.

[53] H. Berendsen, J. Grigera, T. Straatsma, The missing term in effective pair potentials, J. Phys. Chem. 91 (24) (1987) 6269-6271, http://dx.doi.org/10.1021/j100308a038.

[54] M. Szczerba, A. Kuligiewicz, A. Derkowski, V. Gionis, G.D. Chryssikos, A. Kalinichev, Structure and dynamics of water-smectite interfaces: Hydrogen bonding and the origin of the sharp O-Dw/O-Hw infrared band from molecular simulations, Clays Clay Miner. 64 (4) (2016) 452-471, http://dx.doi.org/10.1346/CCMN.2016.0640409.

[55] A. Thompson, H. Aktulga, R. Berger, D. Bolintineanu, W. Brown, P.S. Crozier, P. in 't Veld, A. Kohlmeyer, S. Moore, T. Nguyen, R. Shan, M. Stevens, J. Tranchida, C. Trott, S. Plimpton, Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun. 271 (2022) 10817, http://dx.doi.org/10.1016/j.cpc.2021.108171.

[56] W.M. Brown, P. Wang, S. Plimpton, A. Tharrington, Implementing molecular dynamics on hybrid high performance computers—short range forces, Comput. Phys. Comm. 182 (2011) 898-911, http://dx.doi.org/10.1016/j.cpc.2010.12.021.

[57] W. Brown, A. Kohlmeyer, S. Plimpton, A. Tharrington, Implementing molecular dynamics on hybrid high performance computers—particle-particle particle-mesh, Comput. Phys. Comm. 183 (2012) 449-459, http://dx.doi.org/10.1016/j.cpc.2011.10.012.

[58] R. Hockney, J. Eastwood, Computer Simulation using Particles, Adam Hilger, NY, US, 1988.

[59] W. Swope, H. Andersen, P. Berens, K. Wilson, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, J. Chem. Phys. 76 (1982) 637-649, http://dx.doi.org/10.1063/1.442716.

[60] W. Hoover, Constant-pressure equations of motion, Phys. Rev. A 34 (1986) 2499-2500, http://dx.doi.org/10.1103/PhysRevA.34.2499.

[61] M. Parrinello, A. Rahman, Crystal structure and pair potentials: A molecular-dynamics study, Phys. Rev. Lett. 45 (1980) 1196-1199, http://dx.doi.org/10.1103/PhysRevLett.45.1196.

[62] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method, J. Appl. Phys. 52 (12) (1981) 7182-7190, http://dx.doi.org/10.1063/1.328693.

[63] G.J. Martyna, D.J. Tobias, M.L. Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys. 101 (5) (1994) 4177-4189, http://dx.doi.org/10.1063/1.467468.

[64] W. Shinoda, M. Shiga, M. Mikami, Rapid estimation of elastic constants by molecular dynamics simulation under constant stress, Phys. Rev. B 69 (2004) http://dx.doi.org/10.1103/PhysRevB.69.134103.

[65] J. Cao, G.A. Voth, The formulation of quantum statistical mechanics path centroid density. II. Dynamical properties, J. Chem. Phys. 100 (1993) (1994) 5106-5117, http://dx.doi.org/10.1063/1.467176.

[66] M.E. Tuckerman, B.J. Berne, G.J. Martyna, M.L. Klein, Efficient molecular dynamics and hybrid Monte Carlo algorithms for path integrals, J. Chem. Phys. 99 (August) (1993) 2796-2808, http://dx.doi.org/10.1063/1.465188.

[67] T.D. Kühne, M. Iannuzzi, M. Del Ben, V.V. Rybkin, P. Seewald, F. Stein, T. Laino, R.Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M.H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A.S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G.K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C.J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, J. Hutter, CP2K: An electronic structure and molecular dynamics software package quickstep: Efficient and accurate electronic structure calculations, J. Chem. Phys. 152 (2020) 194103, http://dx.doi.org/10.1063/5.0007045.

[68] S. Göedecker, M. Teter, J. Hütter, Separable dual-space gaussian pseudopotentials, Phys. Rev. B 3 (2) (1996) 1703-1710, http://dx.doi.org/10.1103/PhysRevB.54.1703.

[69] J. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (1996) 3865-3868, http://dx.doi.org/10.1103/PhysRevLett.77.3865.

[70] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132 (2010) 154104, http://dx.doi.org/10.1063/1.3382344.

[71] H. Schlegel, Optimization of equilibrium geometries and transition structures, J. Comput. Chem. 3 (2) (1982) 214-218, http://dx.doi.org/10.1002/jcc.540030212.

[72] D. Rapaport, Hydrogen bonds in water, Mol. Phys. 50 (5) (1983) 1151-1162, http://dx.doi.org/10.1080/00268978300102931.

[73] H. Manzano, A. Ayuela, A. Telesca, P.J.M. Monteiro, J.S. Dolado, Ettringite strengthening at high pressures induced by the densification of the hydrogen bond network, J. Phys. Chem. C 116 (30) (2012) 16138-16143, http://dx.doi.org/10.1021/jp301822e.

[74] N. Stolte, J. Daru, H. Forbert, J. Behler, D. Marx, Nuclear quantum effects in liquid water are marginal for its average structure but significant for dynamics, J. Phys. Chem. Lett. 15 (49) (2024) 12144-12150, http://dx.doi.org/10.1021/acs.jpclett.4c02925.

[75] R. Hill, The elastic behaviour of a crystalline aggregate, Proc. Phys. Soc. Sect. A 65 (5) (1952) 349-354, http://dx.doi.org/10.1088/0370-1298/65/5/307.

[76] A. Togo, L. Chaput, T. Tadano, I. Tanaka, Implementation strategies in phonony and phono3py, J. Phys.: Condens. Matter. 35 (2023) 353001, http://dx.doi.org/10.1088/1361-648X/acd831.

[77] S.M. Mutisya, J.M. de Almeida, C.R. Miranda, Molecular simulations of cement based materials: A comparison between first principles and classical force field calculations, Comput. Mater. Sci. 138 (2017) 392-402, http://dx.doi.org/10.1016/j.commatsci.2017.07.009.

[78] T. Matschei, F. Glasser, Thermal stability of thaumasite, Mater. Struct. 48 (2015) 2277-2289, http://dx.doi.org/10.1617/s11527-014-0309-4.

[79] B. Lothenbach, D.A. Kulik, T. Matschei, M. Balonis, L. Baquerizo, B. Dilnesa, G.D. Miron, R.J. Myers, Cemdata18: A chemical thermodynamic database for hydrated portland cements and alkali-activated materials, Cem. Concr. Res. 115 (2019) 472-506, http://dx.doi.org/10.1016/j.cemconres.2018.04.018.

[80] A. Kalinichev, R. Kirkpatrick, Molecular dynamics modeling of chloride binding to the surfaces of calcium hydroxide, hydrated calcium aluminate, and calcium silicate phases, Chem. Mater. 14 (2002) 3539-3549, http://dx.doi.org/10.1021/cm0107070.

[81] S. Hajilar, B. Shafei, Atomic-scale investigation of physical adsorption of water molecules and aggressive ions to ettringite's surfaces, J. Colloid Interface Sci. 513 (2018) 104-116, http://dx.doi.org/10.1016/j.jcis.2017.09.019.

[82] H. Zhao, Y. Yang, S. Song, X. Shu, W. Chen, W. Jin, Q. Ran, Adsorption conformation of comb-shaped polycarboxylate ethers on ettringite (100) surface: An atomic scales simulation, J. Comput. Biophys. Chem. 20 (1) (2021) 99-108, http://dx.doi.org/10.1142/S2737416521500071.

[83] A. Cody, H. Lee, R. Cody, P. Sprya, The effects of chemical environment on the nucleation, growth, and stability of ettringite [Ca₃Al(OH)₆]₂(SO₄)₃·26H₂O, Cem. Concr. Res. 34 (2004) 869-881, http://dx.doi.org/10.1016/j.cemconres.2003.10.023.

[84] F. Colmenero, A.M. Fernández, O. Almendros-Ginestà, T. Missana, Density functional theory study of the crystal structure and infrared spectrum of a synthesized ettringite mineral, Minerals 14 (2004) 824, http://dx.doi.org/10.3390/min1408024.

[85] H. Heinz, R.A. Vaia, B.L. Farmer, R.R. Naik, Accurate simulation of surfaces and interfaces of face-centered cubic metals using 12-6 and 9-6 Lennard-Jones potentials, J. Phys. Chem. C 112 (44) (2008) 17281-17290, http://dx.doi.org/10.1021/jp801931d.

[86] R.K. Mishra, K. Kanhaiya, J.J. Winetrout, R.J. Flatt, H. Heinz, Force field for calcium sulfate minerals to predict structural, hydration, and interfacial properties, Cem. Concr. Res. 139 (2021) 106262, http://dx.doi.org/10.1016/j.cemconres.2020.106262.

[87] L.G. Baquerizo, T. Matschei, K.L. Scrivener, Impact of water activity on the stability of ettringite, Cem. Concr. Res. 79 (2016) 31-44, http://dx.doi.org/10.1016/j.cemconres.2015.07.008.

[88] T. Honorio, Thermal conductivity, heat capacity and thermal expansion of ettringite and metaetringite: Effects of the relative humidity and temperature, Cem. Concr. Res. 159 (2022) 106865, http://dx.doi.org/10.1016/j.cemconres.2022.106865.

[89] H. van Gog, First-principles study of dehydration interfaces between diaspore and corundum, gibbsite and boehmite, and boehmite and γ-Al₂O₃: Energetic stability, interface charge effects, and dehydration defects, Appl. Surf. Sci. 541 (2021) 148501, http://dx.doi.org/10.1016/j.apsusc.2020.148501.

[90] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO—the Open Visualization Tool, Modelling Simul. Mater. Sci. Eng. 18 (2010) 015012, http://dx.doi.org/10.1088/0965-0393/18/1/015012.
[91] C. Hall, P. Barnes, A.D. Billimore, A. Jupe, X. Turrillas, Thermal decomposition of ettringite Ca₆[Al(OH)₆]₂(SO₄)₃·26H₂O), J. Chem. Soc. Faraday Trans. 92 (12) (1996) 2125–2129, http://dx.doi.org/10.1039/FT9969202125.

[92] L.P. Ogorodova, Y.D. Gritsenko, D.A. Kosova, M.F. Vigasina, L.V. Melchakova, D.A. Ksenofontov, S.K. Dedushenko, Physicochemical and thermochemical study of ettringite, Geochem. Int. 59 (2021) 1188–1197, http://dx.doi.org/10.1134/S0016702921100050.

[93] B. Kostova, V. Petkova, V. Kostov-Kytin, Y. Tzvetanova, G. Avdeev, TG/DTG-DSC and high temperature in-situ XRD analysis of natural thaumasite, Thermochim. Acta 697 (2021) 178863, http://dx.doi.org/10.1016/j.tca.2021.178863.

[94] P. Kostenetskiy, R. Chulkevich, V. Kozyrev, HPC resources of the higher school of economics, J. Phys.: Conf. Ser. 1740 (2021) 012050, http://dx.doi.org/10.1088/1742-6596/1740/1/012050.