# Equation of State, Compressibility, and Vibrational Properties of Brucite over Wide Pressure and Temperature Ranges: Atomistic Computer Simulations with the Modified ClayFF Classical Force Field

Evgeny V Tararushkin, Vasily V Pisarev, Andrey G. Kalinichev

► To cite this version:

Evgeny V Tararushkin, Vasily V Pisarev, Andrey G. Kalinichev. Equation of State, Compressibility, and Vibrational Properties of Brucite over Wide Pressure and Temperature Ranges: Atomistic Computer Simulations with the Modified ClayFF Classical Force Field. Minerals, 2023, 13 (3), pp.408. 10.3390/min13030408 . hal-04074730

HAL Id: hal-04074730

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

Submitted on 15 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

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

Article

# Equation of State, Compressibility, and Vibrational Properties of Brucite over Wide Pressure and Temperature Ranges: Atomistic Computer Simulations with the Modified ClayFF Classical Force Field

Evgeny V. Tararushkin<fnref n="1" />, Vasily V. Pisarev<fnref n="1" />,<fnref n="2" /> and Andrey G. Kalinichev<fnref n="3" />,*

**Abstract:** The behavior of brucite over wide ranges of temperatures and pressures is of great interest for fundamental geochemistry and geophysics. Brucite layers and their octahedral $Mg(OH)_6$ structural units constitute an important structural part of layered dense magnesium hydrous silicates (DMHS), which play a major role in mineral equilibria controlling water balance in the subduction zones of the upper mantle. The ClayFF force field was originally developed for atomistic computer simulations of clays and other layered minerals and their hydrated interfaces. The crystallographic parameters of brucite at 25 °C and 1 bar were used, among several others, to develop the original ClayFF parametrization. Its new recent modification, ClayFF-MOH, can more accurately account for the bending of Mg–O–H angles in the brucite structure, and it was used here to test the applicability of this simple classical model over very wide ranges of temperature and pressure well beyond the range of its original implementation (up to 600 °C and 15 GPa). The pressure and temperature dependencies of brucite crystallographic parameters, the compressibility of the crystal lattice, the coefficients of thermal expansion, and the vibrational spectra were calculated in a series of classical molecular dynamics simulations using the ClayFF-MOH model and compared with a diverse set of available experimental data, including X-ray diffractometry, neutron scattering, IR and Raman spectroscopy. These new results demonstrated that ClayFF-MOH, as simple and approximate as it is, can be quite accurate in predicting many mineral properties at subduction zone conditions, which greatly expands the area of its applicability.

**Keywords:** brucite; atomistic simulations; high pressures; high temperatures; structural properties; elastic properties; bulk modulus; vibrational spectra

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

check for
updates

**Citation:** Tararushkin, E.V.; Pisarev, V.V.; Kalinichev, A.G. Equation of State, Compressibility, and Vibrational Properties of Brucite over Wide Pressure and Temperature Ranges: Atomistic Computer Simulations with the Modified ClayFF Classical Force Field. *Minerals* **2023**, *13*, 408. https://doi.org/10.3390/min13030408

Academic Editors: Lidong Dai, Haiying Hu and Jianjun Jiang

Received: 14 February 2023

Revised: 12 March 2023

Accepted: 13 March 2023

Published: 15 March 2023

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

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

## 1. Introduction

Brucite is a typical layered hydroxide mineral, $M(OH)_2$, where M is a divalent metal ($M = Mg, Ca, Fe, Ni, etc$). Most such hydroxides relate to the layered brucite-type structure [1]. At ambient conditions, brucite crystal structure has a space group $P\bar{3}m1$ symmetry with layers consisting of edge-sharing $[MgO_6]$ octahedra. Hydrogen atoms attached to each oxygen form hydroxyl groups located above and below the octahedral layer and oriented along the c-axis. Despite the fact that the local environment of the OH groups appears to be quite unfavorable for hydrogen bonding between the layers, most authors argue that the nature of the interlayer interactions is due precisely to the hydrogen bonding [2].

The thermodynamic and elastic properties of brucite, $Mg(OH)_2$, especially at high pressures and temperatures, are of great importance in mineralogy and geophysics because brucite itself plays a major role in mineral equilibria controlling water balance in the subduction zones of the Earth's mantle, but also because the octahedral $Mg(OH)_6$ structural

<footnote n="1">¹ International Laboratory for Supercomputer Atomistic Modelling and Multi-Scale Analysis, HSE University, 123458 Moscow, Russia</footnote>

<footnote n="2">² Joint Institute for High Temperatures of the Russian Academy of Sciences, 125412 Moscow, Russia</footnote>

<footnote n="3">³ Laboratoire SUBATECH, UMR 6457-IMT Atlantique, Nantes Université, CNRS/IN2P3, 44307 Nantes, France</footnote>

<footnote>* Correspondence: kalinich@subatech.in2p3.fr</footnote>
units and brucite-type layers constitute an important part of other minerals, such as layered
dense magnesium hydrous silicates (DMHS) important in mantle mineralogy [3–6].

Therefore, the properties of brucite have been extensively studied experimentally by
various techniques over wide ranges of temperature and pressure [2,7–9]. One of the most
interesting crystallographic questions is the behavior of hydrogen (deuterium) atoms and
the nature of hydrogen bonds under these conditions [8,10–13].

The structural and vibrational properties of brucite have been investigated in a number
of experimental studies. In studying the structural parameters of brucite, depending on
pressure and temperature, the positions of atoms and the distances between them are of
interest [11,12,14–16] as well as the unit cell parameters (a, c, *c/a*, *V*) [7–9,11,12,14–19]. Such
structural parameters of brucite were obtained using powder X-ray diffraction [7–9,15,16,18],
single-crystal X-ray diffraction [17] and powder neutron diffraction [20,21]. The mechanical
parameters of brucite are also of great interest, such as isothermal compressibility, bulk
modulus, elastic constants, and thermal expansion coefficient [8,9,15–17,19–22]. Vibration
properties of brucite under high pressure and temperature were also investigated in situ
by Raman spectroscopy [2] and infrared spectroscopy [1,9,23–25]. The most attention has
been paid to the study of OH librations, OH stretching, and the lattice vibrations at various
values of pressure and temperature.

In addition to experimental studies, the structure, and vibrational properties of
brucite were also evaluated by first-principles molecular dynamics simulations based
on the density functional theory (DFT) and by classical molecular dynamic (MD) simula-
tions [2, 6, 10, 12, 13, 26–34]. Such simulations allow data to be obtained on the crystal lattice
parameters, interatomic distances (O–H, O–Mg), elastic constants and bulk modulus, the
coefficients of compressibility and thermal expansion of the brucite crystals, as well as their
vibrational spectra. In addition, DFT modeling allows for the quantitative probing of the
proton disorder [12,30] and dispersive interactions [13] in layered materials. For classical
MD simulations of brucite, the ClayFF force field is most frequently used [26–28]. ClayFF
is a nonbonded force field for classical atomistic simulations of a wide range of clay-related
and cement-related materials [31–37]. In addition to modeling bulk mineral phases, ClayFF
is also well suited for MD simulations of the interaction of water and aqueous solutions
with mineral surfaces [26,28,31–37]. This force field was recently upgraded by adding sev-
eral new terms that allow for a more accurate description of the M–O–H angular bending
in hydroxide crystals [33,34,38].

In this study, we used this modified version of ClayFF to test its applicability in atom-
istic simulations of layered hydroxide minerals at high temperatures and high pressures
by simulating the structural and vibrational properties of brucite. The atomistic model of
brucite was isothermally compressed to 19.86 GPa at 300 K, to 15.15 GPa at 473 K, and to
14.56 GPa at 673 K. The simulation results were then carefully compared with all available
experimental data and with the simulation results with the original version of ClayFF force
field [31]. These new results clearly demonstrated that the ClayFF-MOH model can quite
accurately predict many mineral properties under conditions well beyond its original range
of parametrization [31]. This conclusion greatly expands the area of the applicability for
this simple model, in particular for atomistic simulations of minerals and mineral-water
interactions at subduction zone conditions.

**2. Models and Methods**

2.1. *Structural Models*

A brucite crystal model was constructed based on the structural data from powder
neutron diffraction measurements at room temperature and pressure with the lattice param-
eters corresponding to a trigonal, *P*3*m*1, unit cell (a = b = 3.14979 Å; c = 4.7702 Å; α = β = 90°;
γ = 120°) [21]. This unit cell was used as the initial building block for a simulation su-
percell containing 968 units (11 × 11 × 8 along the *a*, *b*, and *c* crystallographic directions,
respectively). As a result, the simulation supercell contained 4840 atoms and had total
dimensions of 34.65 × 34.65 × 38.16 Å³ (see Figure S1 in the Supplementary Materials).
Periodic boundary conditions were then applied in all three dimensions of the supercell to produce a continuous model of the crystal [39].

## 2.2. Force Field Parameters

All interatomic interactions were calculated using the ClayFF force field [1,34], which has already proven to be highly effective in modeling the structures of many oxides, hydroxides, clays, and other layered materials [33–37]. One of the objectives of the present work was to compare the results obtained by the original ClayFF version [1], further referred to as ClayFF-orig, and its recently developed modification, ClayFF-MOH, [31,33], which more accurately accounts for the bending of M–O–H angles of the hydroxide groups in the crystal by introducing an additional harmonic term to the force field expression:

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

where $\theta_{0,\text{MgOH}} = 110^{\circ}$ and $k_{\text{MgOH}} = 6 \text{ kcal} \cdot \text{mol}^{-1} \cdot \text{rad}^{-2}$ were used for brucite [33]. For a more accurate modeling of the vibrational properties of brucite, the Morse potential was used to describe the hydroxide O–H stretching vibrations instead the original harmonic potential [28,40]:

$$E_{\text{OH}} = D_0 \left[ 1 - e^{-\alpha(r - r_0)} \right]^2 \qquad (2)$$

where $D_0 = 132.2491 \text{ kcal} \cdot \text{mol}^{-1}$, $\alpha = 2.1350 \text{ Å}^{-1}$, $r_0 = 0.9572 \text{ Å}$.

However, the primary objective of the present work was to test the temperature and pressure limits of the ClayFF applicability for the modeling of mineral properties by classical atomistic simulations at thermodynamic conditions beyond the range of original ClayFF parametrization [31].

## 2.3. Simulation Details

LAMMPS simulation package [41,42] was used to perform the MD simulations. As the first step, the potential energy of the brucite model supercell was minimized. Then an MD simulation run in NPT-ensemble was used to bring the supercell to equilibrium at the desired temperature and pressure. The Nosé–Hoover algorithm [39] was employed to control the temperature and pressure in the NPT-simulations. The NPT-simulation run lasted for 0.5 ns using the velocity Verlet’s algorithm [39], which then followed by the NVT-simulation run for another 0.5 ns. The NVT-simulation was used to collect the equilibrium dynamic trajectories of all atoms for further statistical analyses. The coordinates and velocities of all atoms were recorded every 100 fs, which is enough to guarantee the statistical stability of data analysis. A timestep of 1.0 fs was always used for numerical integration of the equations of atomic motions.

## 2.4. Simulation Analysis

Structural and vibrational properties were evaluated in this study by time averaging over the entire simulated equilibrium MD trajectories. In particular, the lattice parameters of the bricute unit cell (a, c, c/a, V) and volume compression of brucite, isothermal bulk moduli, power spectra of vibrational density of states were calculated, as well as maps of atomic density and angular distributions for hydroxyls in the octahedral sheets of brucite.

The unit cell parameters of brucite were directly calculated as averages from the equilibrium NPT-simulations. The isothermal volume compression of the unit cell, $V/V_0$, was also calculated for the corresponding pressures from the equilibrium NPT-simulations.

The isothermal bulk modulus and its pressure derivative were calculated by least squares fitting the third-order Birch–Murnaghan equation of state in the Eulerian finite strain parameter [43,44] to the simulated pressure–volume data:

$$P = \frac{3}{2} K_{T0} \left[ \left( \frac{V_0}{V} \right)^{\frac{7}{3}} - \left( \frac{V_0}{V} \right)^{\frac{5}{3}} \right] \cdot \left\{ 1 + \frac{3}{4} (K'_{T0} - 4) \cdot \left[ \left( \frac{V_0}{V} \right)^{\frac{2}{3}} - 1 \right] \right\} \qquad (3)$$
where $K_{T0}$ and $K'_{T0}$ are the isothermal bulk elastic modulus and its pressure derivative, respectively.

The volumetric thermal expansion coefficient of brucite was obtained by fitting the P-V-T simulation data via constraints on the thermoelastic properties at high-pressure [15,45]. The third-order Birch–Murnaghan equation of state was used to calculate the thermal expansion coefficient through the fitting to a liner relationship [15]:

$$ \alpha(T) = \alpha_0 + \alpha_1 T \qquad (4) $$

The amplitudes of atomic vibrations in the brucite crystal were visualized by calculating the time-averaged contour maps of atomic densities for a slice of the crystal starting from the positions of Mg atoms to the positions of H atoms of the corresponding hydroxyl groups. The thickness of this slice was ~3.0 Å. The atomic density of an atom $\alpha$, $\rho_{\alpha}(x,y)$, was evaluated on the square grid and the time-averaged number of atoms $\langle N_{\alpha}(\Delta x \Delta y) \rangle$ in each grid cell normalized by the surface area of the cell:

$$ \rho_{\alpha}(x, y) = \frac{\langle N_{\alpha}(\Delta x \Delta y) \rangle}{\Delta x \Delta y} \qquad (5) $$

For geometrical analyses of the hydroxyl (O–H bond) angular distributions the angles $\theta$ and $\varphi$ were calculated for both versions of the ClayFF force field, as defined in Figure 1.

<figure><img src="image_5.png" /><figcaption>Figure 1. Angle of slope, $\varphi$, of the structural hydroxyls to the (001) plane of the brucite crystal and the Mg–O–H angle, $\theta$. Green atoms—Mg, red atoms—O, white atoms—H.</figcaption></figure>

The power spectra (PS), providing the quantitative information about the frequencies of atomic vibrations in the system, were obtained as Fourier transformations of the velocity autocorrelation function (VACF) of the corresponding atoms (e.g., [27,37,46,47]):

$$ PS(\omega) = \sum_{t=1}^{n_{\text{VACF}}} \text{VACF}(t) \cdot \cos\left(\frac{\pi}{n_{\text{VACF}}} \cdot t \cdot \omega\right) \qquad (6) $$

where $n_{\text{VACF}}$ is the number of trajectory ‘windows’ used in a particular VACF calculation. The normalized VACFs were calculated from the MD-simulated equilibrium NVT trajectories of the atoms in the modeled system [47]:

$$ VACF(t) = \frac{\sum_{f=1}^{n-n_{\text{VACF}}} \sum_{\text{sel. atoms}} \vec{v}_f \cdot \vec{v}_{f+t}}{\sum_{f=1}^{n-n_{\text{VACF}}} \sum_{\text{sel. atoms}} \vec{v}_f \cdot \vec{v}_f} \qquad (7) $$

where $n$ is the total number of the trajectory ‘windows’; $v_f$ and $v_{f+t}$ are velocities at times $f$ and $f+t$, respectively. They were smoothed by applying a filtering function [47]:

$$ VACF(t) = VACF(t) \cdot \exp\left(-\frac{t}{\tau}\right) \qquad (8) $$
The VACF and PS were calculated for all atoms of the crystal and separately for each atom type (magnesium, oxygen, and hydrogen). We also analyzed the *xy* and *zz* components of the VACF tensor to obtain PS in the *xy* direction (parallel to the plane of brucite layers) and *zz* direction (perpendicular to the layering).

**3. Results and Discussion**

**3.1. Compressibility and Lattice Parameters**

First, we checked how the original and modified versions of the ClayFF force field reproduce the elastic properties of brucite. The compressibility of brucite was studied at temperatures 300, 473, and 673 K and compared with the experimental data [8,9]. For the ClayFF-orig model, the calculations were carried out only up to *T* = 473 K, because at *T* = 673 K, the brucite supercell became already unstable. For the ClayFF-MOH model, the supercell remains stable up to *T* = 673 K, but even with that modification the supercell loses stability at *T* = 873 K.

The compressibility was calculated at pressures up to 20 GPa at *T* = 300 K and up to 15 GPa at 473 K and 673 K. The experimental data and the calculated values of compressibility at *T* = 300 K shown in Figure 2. In this case, good agreement with experimental data is observed in the entire pressure range. A slightly larger deviation is observed for the ClayFF-MOH force field for higher pressures of ~20 GPa. The compressibility of brucite at temperatures of 473 K and 673 K is shown in Figure 2. At 473 K, a good agreement with the experimental data is found again, and the ClayFF-MOH force field results in slightly larger deviations from the experiment. However, only the ClayFF-MOH results are available for *T* = 673 K, as the brucite supercell decomposes with the ClayFF-orig simulations at that temperature even at normal pressure.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Compressibility of brucite: (a) *T* = 300 K; (b) *T* = 473 K and *T* = 673 K. Open circles and rhombs belong to calculated values. Experimental data are from Fei and Mao [7], Fukui et al. [8], Ma et al. [9], and Nagai et al. [16].</figcaption></figure>

A comparison of the crystallographic parameters of brucite at high temperatures and high pressures with experimental data [8] show quite a good agreement (see Figure S2 in the Supplementary Materials). The lattice parameters of the unit cell for ambient conditions calculated in this work virtually coincide with the lattice parameters reported by Pouvreau et al. [33] for the ClayFF-orig version, while for the ClayFF-MOH version of the force field our simulations results (*a* = 3.29 Å, *c* = 4.77 Å, *V* = 44.12 Å³) are slightly higher.
The maximum deviation of the unit cell volume from experimental values is 9% at the highest pressures for both versions of the ClayFF force field at temperatures 300 and 473 K. At *T* = 673 K, the maximum deviation of the unit cell volume for the ClayFF-MOH model is 11%.

The agreement of the parameter *a* of the unit cell with experiment [8] is fairly good, with the discrepancy from the experimental value no more than 5%. There is also a good agreement in the uniaxial compression curve along the *a* axis, including the high-temperature calculations.

For ambient conditions, the deviations from the experimental data [8] for the *c* axis are –1% and –2% for the ClayFF-orig and ClayFF-MOH versions, respectively. However, both versions of the force field do not reproduce the *c* axis temperature dependence. At high pressures, the discrepancy reaches up to 5%. This leads to an incorrect trend of the *c/a* ratio dependence on pressure (Figure S2 in the Supplementary Materials).

## 3.2. Bulk Modulus

The values of bulk modulus *K*<sub>*T*<sub>0</sub></sub> and its pressure derivative *K*'<sub>*T*<sub>0</sub></sub> computed from Equation (3) are shown in Figure 3 in comparison with the experimental data of Fukui et al. [8]. Other experimental data at room temperature, e.g., [15,17,18,21] and the data of DFT(GGA) calculations, e.g., [10,12], are not shown in Figure 3, as they are very close to the data of Fukui et al. [8]. The calculated bulk elastic modulus at temperatures of 300 K and 473 K for both versions of the ClayFF force field shows a good agreement with the experimental data [8], and the maximum discrepancy is only ~3% for the ClayFF-MOH model. For a higher temperature of 673 K, the deviation of the bulk modulus from the experiment for the ClayFF-MOH force field is ~11%, which is also not bad.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a) Bulk elastic modulus of brucite; and (b) pressure derivative of the bulk modulus. Experimental data of Fukui et al. [8] are shown for comparison.</figcaption></figure>

The pressure derivatives of the bulk modulus at *T* = 300 K gave the largest deviation from the experimental data, with the maximum difference of ~21% for the ClayFF-orig model. In fact, the calculated pressure derivatives at room temperature closely agree with the pressure derivatives measured by Catty et al. [21], the largest difference being ~6% for the ClayFF-orig model. With increasing temperature, the first derivative gives a better agreement with experiment [8], and for a *T* = 673 K the difference is ~2% for the ClayFF-MOH model.
### 3.3. Thermal Expansion Coefficient

The values of the two coefficients of thermal expansion ($α_0$ and $α_1$ in Equation (4)) were only calculated from the MD results obtained for the ClayFF-MOH model due to the instability of the brucite supercell at $T = 673$ K for the ClayFF-orig model. The values of the coefficients for the ClayFF-MOH model are $α_0 = 9.04 \times 10^{-5}$ K$^{-1}$ and $α_1 = 13.00 \times 10^{-8}$ K$^{-2}$. The value of $α_0$ is close to the experimentally determined value of $α_0 = 7.3(3) \times 10^{-5}$ K$^{-1}$ [15], while the value of $α_1$ differs more from the experiment ($α_1 = 3.6(13) \times 10^{-8}$ K$^{-2}$ [15]).

### 3.4. Structural Parameters

Nearest-neighbor Mg–O distances were calculated for both ClayFF versions (Figure S3a). At ambient conditions, both versions of the force field give slightly larger Mg–O distances (2.2–2.3 Å) compared to the experimental [21] and DFT [28] results (2.12–2.13 Å). With increasing pressure, at $T = 300$ K, the Mg–O distance reduces, as well as its mean-squared deviation. Overall, the agreement in the Mg–O distances between the ClayFF models and experimental data [16,21] even improves at higher pressures. The same trends of decreasing Mg–O distance and its fluctuations with increasing pressure are maintained at higher temperatures (473 and 673 K). The distances obtained by both ClayFF models virtually coincide, the difference is only in a larger fluctuation for the ClayFF-orig version, which indicates a more stable brucite structure when the M–O–H angles in the ClayFF-MOH model are more accurately constrained.

The distributions and averages of the hydroxyl O–H distances in brucite for both versions of the ClayFF force field were also calculated. At ambient pressure, the O–H distance is close to the O–H bond length corrected for thermal motion [21] and the bond lengths calculated by DFT [10]. The average O–H distance slightly increases with increasing pressure, which does not agree with experimental data [21] but agrees well with the DFT results at high pressures [10]. The tendency to increase the O–H distance with increasing pressure persists at elevated temperatures. Similar to the Mg–O distance, increasing temperature leads to a slight increase of the average O–H bond length and its dispersion. The lack of good agreement with these experimental data at elevated pressures can be explained by the behavior of the Morse potential used here for modeling the O–H bond of hydroxyls, which has a fairly strong repulsive branch of the potential when hydrogen atom approaches the oxygen. On the other hand, an increase in O–D distance with increasing pressures is indeed experimentally observed for deuterated brucite [20]. So, a better judgement about the quality of the description of this parameter with the ClayFF model could probably be made when more experimental data are available. All the comparisons are shown in Figure S3b in the Supplementary Materials.

The angular distribution of the hydroxyl orientations relative to the brucite (001) plane of layering at $T = 300$ K and high pressures for both versions of the ClayFF model are shown in Figure 4. These distributions are qualitatively different between the two ClayFF versions. At ambient pressure, the direction normal to the brucite layers is the most probable hydroxyl orientation. However, the distribution is significantly broader for the ClayFF-orig model where hydroxyl orientations are not constrained and have the ability to lean toward the brucite layer. With increasing pressure, the orientation normal to the layering remains the most probable for the ClayFF-MOH model, but the earlier ClayFF-orig model demonstrates the appearance of another probability maximum. For example, the most probable orientation is ~63° at $P = 10.35$ GPa. Small deviations of hydroxyl orientations from the normal to the brucite surfaces at ambient pressure is in good agreement with the angles calculated from experimental data [48]. According to experimental data [21], angle $φ$ (see Figure 1) deviates from the normal by only 6.8° at $P = 10.9$ GPa, compared to the same angle under ambient conditions, which, in turn, indicates that hydroxyls remain nearly normal to the layering even at very high pressures, so that the angular distributions for the ClayFF-orig model exhibit unphysical behavior under these conditions. At higher temperatures, in all cases, the overall shape of the hydroxyl orientational distribution remains the same (Figure 4b), but a broadening
is observed, consistent with higher thermal energies of angular vibrations at higher temperatures.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** (a) Angular distributions of hydroxyl orientations at *T* = 300 K and high pressures; and (b) angular distributions of hydroxyl orientations at *T* = 473 K and *T* = 673 K at *P* = 1 bar.</figcaption></figure>

The nature of the second maximum at ~60° for the ClayFF-orig model can be understood by considering the distribution of atomic densities projected onto the brucite *ab* crystallographic plane. Atomic density maps for brucite layer at *T* = 473 K and *P* = 15.15 GPa for both ClayFF versions are shown in Figure 5.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Contour maps of atomic densities at *T* = 473 K and *P* = 15.15 GPa for 3 × 3 unit cell. (a) ClayFF-orig; and (b) ClayFF-MOH. Green contours—Mg, red contours—O, black contours—H.</figcaption></figure>

The ClayFF-orig model shows three predominant orientations of the O–H groups forming an equilateral triangle of probability distributions. These orientations are not accidental, they emerge because the H atoms of the unconstrained O–H groups can easily form hydrogen bonds with one of the three O atoms of the opposite brucite layer. For the ClayFF-MOH model, the atomic density maps look quite different, the O–H groups are now constrained in their orientations and their H atoms are located more or less above the respective O atoms. In this case, hydrogen bonds with the opposite layer are also formed, but instead of one strong bond towards one of the opposing O atoms, three weaker bonds with three nearest O atoms exist simultaneously. In other words, due to the energy penalty from the new M–O–H angle bending terms of the force field, this interaction is not chaotic and H atom stays near the center of a triangle formed by the closest opposite-layer O atoms.
### 3.5. Power Spectra

The spectra of vibrational density of states, or so-called power spectra (Equations (6)–(8)) were calculated for both ClayFF versions and shown in Figure 6. Experimental studies [1,2,9,14,23–25] and atomistic simulations [6,10,13,27,28] indicate that the most intense vibrational bands of brucite under ambient conditions are in the frequency ranges of 270–900 cm⁻¹ and 3550–3890 cm⁻¹. The lower frequency bands correspond to the O–H librational modes (E_g and E_u) and to the lattice vibrations (A₁g, A₂u, E_g and E_u). The higher frequency bands correspond to the O–H stretching modes (A₁g and A₂u) [1,2,9,10,13,14,25,27].

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Power spectra of brucite at T = 298 K and P = 1 bar: (a) full spectrum; and (b) lower frequency modes.</figcaption></figure>

The calculated power spectra of brucite under ambient conditions are in good agreement with the previous results [28] for the ClayFF force field with the constrained Mg–O–H angles. Zeitler et al. [28] have already conducted a detailed analysis of the peak positions in the low-frequency range of the spectrum. In particular, they have shown that taking into account the Mg–O–H angle allows for a slightly better agreement between the calculated classical MD and DFT results for the O–H librational modes, as it was also later confirmed by Pouvreau et al. [33]. Our decomposition of the total vibrational spectra of brucite into the contributions due to O and H atomic motions in the *xy* plane also shows the difference in the O–H librational modes between the two ClayFF versions (Figure S4 in the Supplementary Materials).

Our calculations show the following features of the spectral peaks in the lower frequency range at high pressures and temperatures:

* most peaks exhibit a red shift, and their intensity decreases with increasing temperature (Figure S5 in the Supplementary Materials);

* the peaks exhibit blue shift, and their intensity decreases with increasing pressure (Figure S6 in the Supplementary Materials).

Such tendencies at high pressures agree with the experimental data [2]. Also, at high pressures there is a merging and splitting of some peaks (Figure S6 in the Supplementary Materials).

The higher frequency range of the spectra describes the O–H bond-stretching modes, which correspond to the modes A₁g (Raman) and A₂u (infrared) in experiments. For these modes, the peaks are located at 3573 cm⁻¹ and 3612 cm⁻¹ for the ClayFF-orig and ClayFF-MOH models, respectively (Figure 6). As the frequency is mainly determined by the O–H bond parameters, the peak positions are in good agreement with experimental studies and
DFT calculations [1,2,6,9,10,13–15,24,27] for both ClayFF versions. However, the ClayFF versions differ in the width and intensity of the peaks: the O–H stretching peak is narrower and more intense for the ClayFF-MOH, compared to ClayFF-orig model.

At ambient pressure, the O–H stretching peak does not change its position with tem-
perature in the ClayFF-MOH model. The increasing temperature broadens the peak and
reduces its intensity due to the thermal motion (Figure S5 in the Supplementary Materials).
However, the peaks shift with the increasing pressure indicating a change in the effective
bond stiffness due to the interaction with the neighboring atoms. The computed O–H
stretching peak positions for both versions of the ClayFF force field are shown in Figure 7
together with Raman experimental data and DFT results. One can observe that the change
of the peak position is smaller for ClayFF-MOH than for ClayFF-orig. Such a difference is
likely caused by the difference in the orientational dynamics of the O–H bonds between
the two force field versions: as discussed in Section 3.4, the O–H bonds in ClayFF-orig
are inclined towards a specific oxygen atom, while in ClayFF-MOH they stay more or
less normal to the surface. That difference in the local atomic environments produces the
difference in the stiffness of the O–H groups, in addition to the bond potential. Overall,
ClayFF-MOH shows a better agreement with the experimental data [2] and DFT calcula-
tions [10], especially at the highest pressures (~15 GPa and ~20 GPa), showing that the
constrained O–H orientation leads to a better description of the vibrational dynamics. The
above-described behavior of the peaks at high pressures persists at elevated temperatures
as well (Figure S7 in the Supplementary Materials).

<figure><img src="image_11.png" /><figcaption>Figure 7. Calculated and experimental O–H stretching peak positions at high pressures and room temperature. Experimental data [2] and DFT results [10] are shown for comparison.</figcaption></figure>

**4. Conclusions**

The properties of brucite under ambient conditions and up to high temperature and
high-pressure conditions were investigated by classical MD simulations using the original
ClayFF force field [31] and its more recent modification ClayFF-MOH [34], where the
bending motions of the Mg–O–H angles are more accurately constrained by additional
terms based on the DFT calculations [33]. For both models, the crystallographic, elastic,
and vibrational properties acceptably agreed with available experimental data and the
results of more rigorous DFT calculations. However, the introduction of the Mg–O–H angle
bending term made it possible to simulate brucite at higher temperatures and pressures
(T = 673 K and P = 15 GPa) than with the original ClayFF parameterization. The isothermal
compressibility of brucite was also reproduced well by both versions of the force field,
including at high temperatures. The bulk modulus was also well reproduced, while at
higher temperatures the agreement with the experimental data worsened. The crystal
lattice unit cell parameters (a, c, V and c/a) were reproduced with acceptable accuracy. In
most cases (a, c and V), the trends of the lattice parameters' dependence on pressure and
temperature for the ClayFF model agreed well with the experimental data and DFT results.
The application of the new more accurate ClayFF-MOH force field results in the formation of a more realistic hydrogen bonding structure between brucite layers, which leads to a more stable brucite structure at higher temperatures and pressures, and a more realistic description of its vibrational spectrum.

It is important to emphasize that both versions of the ClayFF model were intentionally used here well beyond the range of temperatures and pressures for which this model was initially parameterized [31,33,34]. Therefore, the results of this work are quite encouraging in view of the applicability of these simple models to the classical atomistic simulations of more complex hydrous minerals under the subduction zone conditions (e.g., [49,50]). This will be the focus of our further investigation.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/min13030408/s1. Figure S1: Two projections of the brucite simulation supercell; Figure S2: (a) Volume of the unit cell for brucite as function of pressure; (b) pressure dependence of the parameter *a* of brucite; (c) pressure dependence of the parameter *c* for brucite; (d) ratio of unit cell parameters *c*/*a* for brucite as a function of pressure; Figure S3: (a) Pressure dependence of the simulated Mg–O distances in [MgO₆] octahedra; (b) pressure dependence of the simulated O–H distances in structural hydroxyls; Figure S4: Power spectra for Hₓystable for both versions of the ClayFF force field; Figure S5: Power spectra of brucite at *P* = 1 bar and high temperatures for the ClayFF-MOH model; Figure S6: Power spectra of brucite at *T* = 300 K and high pressures for the ClayFF-MOH model; Figure S7: Power spectra of brucite at *T* = 673 K and high pressures for the ClayFF-MOH model.

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

**Funding:** This study was sponsored by Agence Nationale pour la Gestion des Déchets Radioactifs and National Research University Higher School of Economics.

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

**Acknowledgments:** This research was funded by the HSE University Basic Research Program and was supported by computational resources of HPC facilities at HSE University [51]. 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. The APC was funded by publication vouchers from previous MDPI reviewing activity.

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

**References**

1. Lutz, H.D.; Moller, H.; Schmidt, M. Lattice vibration spectra. Part LXXXII. Brucite-type hydroxides M(OH)₂ (M = Ca, Mn, Co, Fe, Cd)—IR and Raman spectra, neutron diffraction of Fe(OH)₂. *J. Mol. Struct.* **1994**, *328*, 121–132. [CrossRef]

2. Duffy, T.S.; Meade, C.; Fei, Y.; Mao, H.-K.; Hemley, R.J. High-pressure phase transition in brucite, Mg(OH)₂. *Am. Mineral.* **1995**, *80*, 222–230. [CrossRef]

3. Thompson, A.B. Water in the Earth’s upper mantle. *Nature* **1992**, *358*, 295–302. [CrossRef]

4. Kawamoto, T. Hydrous phases and water transport in the subducting slab. *Rev. Mineral. Geochem.* **2006**, *62*, 273–289. [CrossRef]

5. Ivanov, A.V.; Litasov, K.D. The deep water cycle and flood basalt volcanism. *Int. Geol. Rev.* **2013**, *56*, 1–14. [CrossRef]

6. Hermann, A.; Mookherjee, M. High-pressure phase of brucite stable at Earth’s mantle transition zone and lower mantle conditions. *Proc. Nat. Acad. Sci. USA* **2016**, *113*, 13971–13976. [CrossRef]

7. Fei, Y.; Mao, H.-K. Static compression of Mg(OH)₂ to 78 GPa at high temperature and constraints on the equation of state of fluid H₂O. *J. Geophys. Res. Solid Earth* **1993**, *98*, 11875–11884. [CrossRef]

8. Fukui, H.; Ohtaka, O.; Suzuki, T.; Funakoshi, K. Thermal expansion of Mg(OH)₂ brucite under high pressure and pressure dependence of entropy. *Phys. Chem. Miner.* **2003**, *30*, 511–516. [CrossRef]

9. Ma, M.; Liu, W.; Chen, Z.; Liu, Z.; Li, B. Compression and structure of brucite to 31 GPa from synchrotron X-ray diffraction and infrared spectroscopy studies. *Am. Mineral.* **2013**, *98*, 33–40. [CrossRef]

10. Mitev, P.D.; Gajewski, G.; Hermansson, K. Anharmonic OH vibrations in brucite: Small pressure-induced redshift in the range 0–22 GPa. *Am. Mineral.* **2009**, *94*, 1687–1697. [CrossRef]
11. Xu, H.; Zhao, Y.; Hickmott, D.D.; Lane, N.J.; Vogel, S.C.; Zhang, J.; Daemen, L.L. High-temperature neutron diffraction study of deuterated brucite. *Phys. Chem. Miner.* **2013**, *40*, 799–810. [CrossRef]

12. Mookherjee, M.; Stixrude, L. High-pressure proton disorder in brucite. *Am. Mineral.* **2006**, *91*, 127–134. [CrossRef]

13. Ugliengo, P.; Zicovich-Wilson, C.M.; Tosoni, S.; Civalleri, B. Role of dispersive interactions in layered materials: A periodic B3LYP and B3LYP-D* study of Mg(OH)₂, Ca(OH)₂ and kaolinite. *J. Mater. Chem.* **2009**, *19*, 2564–2572. [CrossRef]

14. Chakoumakos, B.C.; Loong, K.K.; Schultz, A.J. Low-temperature structure and dynamics of brucite. *J. Phys. Chem. B* **1997**, *101*, 9458–9462. [CrossRef]

15. Xia, X.; Weidner, D.J.; Zhao, H. Equation of state of brucite; single-crystal Brillouin spectroscopy study and polycrystalline pressure-volume-temperature measurement. *Am. Mineral.* **1998**, *83*, 68–74. [CrossRef]

16. Nagai, T.; Ito, T.; Hattori, T.; Yamanaka, T. Compression mechanism and amorphization of portlandite, Ca(OH)₂: Structural refinement under pressure. *Phys. Chem. Miner.* **2000**, *27*, 462–466. [CrossRef]

17. Duffy, T.S.; Shu, J.; Mao, H.-K.; Hemley, R.J. Single-crystal X-ray diffraction of brucite to 14 GPa. *Phys. Chem. Miner.* **1995**, *22*, 277–281. [CrossRef]

18. Horita, J.; dos Santos, A.M.; Tulk, C.A.; Chakoumakos, B.C.; Polyakov, V.B. High-pressure neutron diffraction study on H-D isotope effects in brucite. *Phys. Chem. Miner.* **2010**, *37*, 741–749. [CrossRef]

19. Redfern, S.A.T.; Wood, B.J. Thermal expansion of brucite, Mg(OH)₂. *Am. Mineral.* **1992**, *77*, 1129–1132.

20. Parise, J.B.; Leinenweber, K.; Weidner, D.J.; Tan, K.; Von Dreele, R.B. Pressure-induced H bonding: Neutron diffraction study of brucite, Mg(OD)₂, to 9.3 GPa. *Am. Mineral.* **1994**, *79*, 193–196.

21. Catti, M.; Ferraris, G.; Hull, S.; Pavese, A. Static compression and H disorder in brucite, Mg(OH)₂, to 11 GPa: A powder neutron diffraction study. *Phys. Chem. Miner.* **1995**, *22*, 200–206. [CrossRef]

22. Jiang, F.; Speziale, S.; Duffy, T.S. Single-crystal elasticity of brucite, Mg(OH)₂, to 15 GPa by Brillouin scattering. *Am. Mineral.* **2006**, *91*, 1893–1900. [CrossRef]

23. Kruger, M.B.; Williams, Q.; Jeanloz, R. Vibrational spectra of Mg(OH)₂ and Ca(OH)₂ under pressure. *J. Chem. Phys.* **1989**, *91*, 5910–5915. [CrossRef]

24. de Oliveira, E.F.; Hase, Y. Infrared study and isotopic effect of magnesium hydroxide. *Vib. Spectrosc.* **2001**, *25*, 53–56. [CrossRef]

25. Zhu, X.; Guo, X.; Smyth, J.R.; Ye, Y.; Wang, X.; Liu, D. High-temperature vibrational spectra between Mg(OH)₂ and Mg(OD)₂: Anharmonic contribution to thermodynamics and D/H fractionation for brucite. *J. Geophys. Res. Solid Earth* **2019**, *124*, 8267–8280. [CrossRef]

26. Wang, J.W.; Kalinichev, A.G.; Kirkpatrick, R.J. Molecular modeling of water structure in nano-pores between brucite (001) surfaces. *Geochim. Cosmochim. Acta* **2004**, *68*, 3351–3365. [CrossRef]

27. Braterman, P.S.; Cygan, R.T. Vibrational spectroscopy of brucite: A molecular simulation investigation. *Am. Mineral.* **2006**, *91*, 1188–1196. [CrossRef]

28. Zeitler, T.R.; Greathouse, J.A.; Gale, J.D.; Cygan, R.T. Vibrational analysis of brucite surfaces and the development of an improved force field for molecular simulation of interfaces. *J. Phys. Chem. C* **2014**, *118*, 7946–7953. [CrossRef]

29. Ulian, G.; Valdre, G. Equation of state and second-order elastic constants of portlandite Ca(OH)₂ and brucite Mg(OH)₂. *Phys. Chem. Miner.* **2018**, *46*, 101–117. [CrossRef]

30. Schaack, S.; Depondt, P.; Huppert, S.; Finocchi, F. Quantum driven proton diffusion in brucite-like minerals under high pressure. *Sci. Rep.* **2020**, *10*, 8123. [CrossRef]

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

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

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

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

35. Kalinichev, A.G.; Wang, J.W.; Kirkpatrick, R.J. Molecular dynamics modeling of the structure, dynamics and energetics of mineral-water interfaces: Application to cement materials. *Cem. Concr. Res.* **2007**, *37*, 337–347. [CrossRef]

36. Mishra, R.K.; Mohamed, A.K.; Geissbühler, D.; Manzano, H.; Jamil, T.; Shahsavari, R.; Kalinichev, A.G.; Galmarini, S.; Tao, L.; Heinz, H.; et al. cemff: A force field database for cementitious materials including validations, applications and opportunities. *Cem. Concr. Res.* **2017**, *102*, 68–89. [CrossRef]

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

38. Pouvreau, M.; Greathouse, J.A.; Cygan, R.T.; Kalinichev, A.G. Structure of hydrated kaolinite edge surfaces: DFT results and further development of the ClayFF classical force field with Metal-O-H angle bending terms. *J. Phys. Chem. C* **2019**, *123*, 11628–11638. [CrossRef]

39. Allen, M.P.; Tildesley, D.J. *Computer Simulation of Liquids*, 2nd ed.; Oxford University Press: New York, NY, USA, 2017. [CrossRef]
40. Greathouse, J.A.; Durkin, J.S.; Larentzos, J.P.; Cygan, R.T. Implementation of a Morse potential to model hydroxyl behaviour in phyllosilicates. *J. Chem. Phys.* **2009**, *130*, 134713. [CrossRef]

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

42. Thompson, A.P.; Aktulga, H.M.; Berger, R.; Bolintineanu, D.S.; Brown, W.M.; Crozier, P.S.; in't Veld, P.J.; Kohlmeyer, A.; Moore, S.G.; Nguyen, T.D.; et al. LAMMPS—a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comp. Phys. Comm.* **2022**, *271*, 108171. [CrossRef]

43. Birch, F. Finite strain isotherm and velocities for single-crystal and polycrystalline NaCl at high pressures and 300°K. *J. Geophys. Res. Solid Earth* **1978**, *83*, 1257–1268. [CrossRef]

44. Katsura, T.; Tange, Y. A simple derivation of the Birch-Murnaghan equations of state (EOSs) and comparison with EOSs derived from other definitions of finite strain. *Minerals* **2019**, *9*, 745. [CrossRef]

45. Jackson, I.; Rigden, S.M. Analysis of P-V-T data: Constraints on the thermoelastic properties of high-pressure minerals. *Phys. Earth Planet. Inter.* **1996**, *96*, 85–112. [CrossRef]

46. Kirkpatrick, R.J.; Kalinichev, A.G.; Wang, J.; Hou, X.; Amonette, J.E. Molecular modeling of the vibrational spectra of interlayer and surface species of layered double hydroxides. *CMS Workshop Lect.* **2005**, *13*, 239–285. [CrossRef]

47. Szczerba, M.; Kuligiewicz, A.; Derkowski, A.; Gionis, V.; Chryssikos, G.D.; Kalinichev, A.G. 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 Min.* **2016**, *64*, 452–471. [CrossRef]

48. Desgranges, L.; Calvarin, G.; Chevrier, G. Interlayer interactions in M(OH)₂: A neutron diffraction study of Mg(OH)₂. *Acta Crystallogr. Sect. B* **1996**, *52*, 82–86. [CrossRef]

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

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

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

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