Investigating the effect of Cu2+ sorption in
montmorillonite using density functional theory and
molecular dynamics simulations

Yalda Pedram, Yaoting Zhang, Scott Briggs, Chang Seok Kim, Laurent

Brochard, Andrey G Kalinichev, Laurent Karim Béland

► To cite this version:

Yalda Pedram, Yaoting Zhang, Scott Briggs, Chang Seok Kim, Laurent Brochard, et al. Investigating the effect of Cu2+ sorption in montmorillonite using density functional theory and molecular dynamics simulations. Comput.Mater.Sci., 2025, 248, pp.113554. 10.1016/j.commatsci.2024.113554. hal-04888148

HAL Id: hal-04888148

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

Submitted on 6 Jan 2026

**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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Full length article

# Investigating the effect of Cu²⁺ sorption in montmorillonite using density functional theory and molecular dynamics simulations

Yalda Pedram<sup>a ID</sup>, Yaoting Zhang<sup>a ID</sup>, Scott Briggs<sup>b ID</sup>, Chang Seok Kim<sup>b ID</sup>, Laurent Brochard<sup>c ID</sup>, Andrey G. Kalinichev<sup>d ID</sup>, Laurent Karim Béland<sup>a ID,*</sup>

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

## ARTICLE INFO

### Keywords:

Montmorillonite

Density functional theory

Molecular dynamics

ClayFF

Swelling pressure

Cation mobility

## ABSTRACT

Montmorillonite (MMT) is the main mineral component of bentonite, which is currently proposed as a sealing material in deep geological repositories (DGRs) for used nuclear fuel. In the Canadian program, which will utilize copper-cladded used fuel containers, safety analysis considers the effect of copper corrosion, during which Cu²⁺ ions could potentially be adsorbed by the surrounding MMT. In such a scenario, ion exchange between Na⁺ and Cu²⁺ is expected. In this study, a multiscale approach that combines electronic density functional theory (DFT) and force-field-based molecular dynamics (MD) simulations was employed to study the effect of introducing Cu²⁺ ions to MMT. An extension to the ClayFF force field is parameterized and validated using DFT to model how Cu²⁺ interacts with clay systems. MD simulations were performed to calculate the interaction free energies between MMT platelets containing Cu²⁺ ions (Cu-MMT) and compared them to inter-platelet interaction energies in Na-MMT and Ca-MMT. Our calculations suggest Cu-MMT develops swelling pressures between those of Ca-MMT and Na-MMT. Furthermore, our MD simulations suggest that Cu²⁺ has MMT interlayer mobility that is significantly slower than that of Ca²⁺.

## 1. Introduction

Bentonite is used by diverse industries [1,2] because of its high swelling potential, strong adsorption capacity, temperature durability, and low permeability [3–5]. It is added to drilling fluids within the oil and gas exploration sector [6–8], acts as an impermeable barrier for geotechnical engineering projects to prevent water infiltration [2,9], and is used in the ceramics industry as a binding agent [2]. In the pharmaceutical and medical domains, it is utilized as a drug delivery systems, as a binder in tablet formulations and for absorbing poisons [2, 10, 11]. In environmental engineering, bentonite is used as efficient adsorbent of heavy metals and contaminants during wastewater treatment processes [2,12,13].

Bentonite is also under consideration in the context of nuclear waste management, being proposed as a sealing material in deep geological repositories (DGRs) [2,14,15]. Proposed DGRs involve multiple barriers, typically including a radioactive waste form enclosed within a long-lived container, surrounded by bentonite-based materials. In the Canadian design, the copper-coated canisters within this system will undergo changing environmental conditions, transitioning from an oxic

environment over a few weeks or months to anoxic conditions thereafter [16,17]. During the initial repository stage, oxygen will interact with copper, resulting in a surface film of Cu₂O or CuO/Cu(OH)₂ [18]. Depending on groundwater chemistry, it is possible that some Cu ions would dissolve in the surrounding solution, likely in the form of Cu²⁺ or possibly as an anionic complex such as CuCl₂⁻ in high chloride systems and would interact with the surrounding bentonite. This investigation focuses on the cationic species Cu²⁺.

The primary component of bentonite is the montmorillonite (MMT) mineral. MMT consists of silicon tetrahedral (T) and aluminum/magnesium octahedral layers (O), forming TO and TOT layers at the nanoscale, with hydrated interlayers and hydroxyl group edges [19,20]. A distinguishing feature of MMT is its permanent net negative charge, resulting from isomorphic metal substitution. For instance, Al³⁺ or Fe³⁺ can replace Si⁴⁺ in T layers, while Mg²⁺ or Li⁺ can substitute Al³⁺ in O layers [20,21]. Layer-charged MMT is balanced by hydrated cations, such as Ca²⁺ and Na⁺. The swelling mechanism of these minerals is driven by the hydration of adsorbed ions between the layers due to their

<footnote><sup>a</sup> Department of Mechanical and Materials Engineering, Queen's University, Kingston, ON K7L 3N6, Canada</footnote>

<footnote><sup>b</sup> Nuclear Waste Management Organization, Toronto M4T 2S3, Ontario, Canada</footnote>

<footnote><sup>c</sup> Laboratoire Navier, Ecole des Ponts ParisTech, Université Gustave Eiffel, CNRS, Marne-la-Vallée, France</footnote>

<footnote><sup>d</sup> Laboratoire SUBATECH, UMR 6457–Institut Mines Télécom Atlantique, Nantes Université, CNRS/IN2P3, 44307 Nantes, France</footnote>

<footnote>* Corresponding author.
E-mail address: laurent.beland@queensu.ca (L.K. Béland).</footnote>
high hydration energy, leading to charge neutralization; such swelling clays are also known as smectite clays [20,22,23].

Several experiments have shown two fundamental processes governing MMT behavior: crystalline swelling and osmotic swelling [24–26]. Crystalline swelling, linked to hydrated interlayer cations, is stepwise and confined within a basal spacing range of ~9–20 Å which is observed in a wide range of clay minerals [22,24,27]. In contrast, osmotic swelling is a continuous process, involving a range greater than ~35–40 Å, and results from ion concentration disparities, also increases linearly with increasing water content [22,24]. In the osmotic swelling range, Na-montmorillonite (Na-MMT) can absorb water, expanding significantly beyond the limits of crystalline swelling [22]. Norrish et al. [25] reported 20-fold volume expansion of Na-MMT in both osmotic and crystalline regions. This expansion is driven by a combination of electrostatic and van der Waals forces, as well as pressure from cations between the clay layers [24,25]. Of note, the extent of swelling in Na-MMT is linked to multiple factors, including dry density, solution salinity, and temperature [28].

There exists a transition region known as the 'Ravina and Low' zone [29], which marks the abrupt shift between these two swelling regimes. However, this transition, referred to as *forbidden* by Ravina and Low [29], has not been observed in X-ray diffraction experiments [25,26,30,31]. Fink et al. [26] demonstrated that this forbidden layer occurs in Wyoming Na-MMT, where isomorphic substitution occurs in both the tetrahedral and octahedral layers, and not in Na-MMTs with only octahedral substitutions. Furthermore, Meleshyn et al. [32] employed Monte Carlo simulations to reveal that interlayer spaces within such Na-MMTs can be locked due to a chainlike structure composed of substituted tetrahedron-Na cation-water molecule-Na cation-substituted tetrahedron' between layers. This space can remain locked until the water content approaches 850 mg/g.

Classical DLVO theory [33], named after its developers Derjaguin, Landau, Verwey, and Overbeek, describes the interactions between colloidal particles in a liquid medium. The theory considers two colloidal particles and characterizes their interaction as a balance between repulsive electrostatic forces and attractive van der Waals forces. However, the DLVO theory fails to capture nanoscale forces between MMT platelets and the influence of divalent cations which makes it inadequate for understanding swelling behavior in the range of 1–2 nm [34–36]. Molecular dynamics (MD) simulations overcome these limitations and allow to better understand MMT at these very fine scales. Atomistic MD simulations directly probe nanoscale interactions and are a good complement to traditional characterization methods such as scanning and transmission electron microscopy and X-ray diffraction (XRD) [37].

ClayFF has been widely used in the literature for the atomistic simulations of clays and other mineral structures. Of particular interest for the present study, MMT is routinely modeled using the ClayFF force-field, developed by Randall T. Cygan and collaborators [38]. This force-field is based on electrostatic and Lennard-Jones pair interactions and will be elaborated upon in the method section. For example, Ho et al. [39] used a modified ClayFF force-field, as discussed later, to explore the impact of hydroxyl edges on Na-MMT hydration. Hydroxyl groups hindered the entry of the initial water layer (1-W) to the interlayer. The study indicated an alignment of the 0-W state at 9.85 Å, consistent with XRD results for dry MMT (9.6–9.8 Å) [40]. The study clarified the role of hydrogen bonds between hydroxyl layers as gatekeepers for water entry between interlayers. Sun et al. [41] used a ClayFF-based model to explore swelling pressure in semi-periodic Na-MMT sheets in various solutions; swelling pressure increased exponentially with dry density, regardless of the solution, with the highest swelling pressure observed in pure water. In another MD-based study, Akinwunmi et al. [42] examined the swelling characteristics of periodic Ca-MMT sheets and periodic Na-MMT sheets in pure water. According to the MD simulations, at low to moderate dry densities, Ca-MMT's swelling pressure was lower than Na-MMT's. At high dry density, Ca-MMT exerted greater swelling pressure than Na-MMT [43].

In order to study the influence of Cu²⁺ on MMT at the atomistic level, existing ClayFF force-field need to be extended to handle Cu²⁺ in clay. To do so, electronic density functional theory (DFT), which is an approximation to quantum mechanics, will be employed as a benchmark to fit this extended force-field. Our main objectives are twofold:

1. Parameterize and validate ClayFF extension for Cu²⁺: We aim to develop and validate an extension of ClayFF with new interaction parameters for Cu²⁺ based on DFT calculations.

2. Investigate the comparative behavior of counter-ions in MMT: We seek to understand the interaction and behavior of MMT platelets with different counter-ions (Na⁺, Ca²⁺, and Cu²⁺).

To achieve these objectives, we performed MD simulations to calculate interaction free energies between MMT platelets containing Cu²⁺ and Ca²⁺ counter-ions and compared these to interaction energies with Na⁺ counter-ions. In addition, we estimated the swelling pressure between MMT platelets and compared the cases of Na⁺, Ca²⁺, and Cu²⁺. Finally, we computed the self-diffusion coefficient of these three cations when confined in the interlayer space.

This work is divided into the following sections: in section Method we will explain the computational procedures, simulation parameters and employed MMT models. Section Results and analysis is divided into three sub-sections. First, we will discuss the results of fitting procedures, then the validation of the extended ClayFF force-field and finally the obtained PMFs, swelling pressures and mobility of cations will be shown. 2.

## 2. Method

ClayFF is not designed to handle Cu²⁺ on MMT. To overcome this limitation we conducted a series of electronic density functional theory (DFT) simulations. We will begin by discussing the DFT setup. Following that, we will discuss the ClayFF force-field, which forms the basis of our simulations. The MD setup utilized to investigate the behavior of our systems will also be discussed, including an introduction to the potential of mean force (PMF) approach employed to calculate interaction free energies.

### 2.1. *Density functional theory*

DFT simulations were performed as implemented in the open-source software Quantum Espresso [44] with projector-augmented waves (PAW) [45–47]. The generalized gradient approximation (GGA) by Perdew–Burke–Ernzerhof (PBE) was employed to calculate the exchange–correlation energy [48]. To account for long-range van der Waals interactions, the D3 dispersion correction method of Grimme was applied [49], which improves the accuracy of interaction energies, particularly in layered materials like MMT. We prepared a periodic Arizona MMT containing 2 × 1 × 1 unit cells with the size of 10.54 ×9.15×6.75 Å using the ATOM package [50] with a total charge of −2. To neutralize the negative charge of the MMT, we introduced one Cu²⁺, resulting in Cu-MMT with the formula [Cu].[Si₁₆[Mg₂Al₆]O₄₀(OH)₈]. A similar formula was used for Ca-MMT by substituting Cu²⁺ with Ca²⁺.

In another scenario, similar structures were used in separate simulations with three water molecules, two Na⁺, one Ca²⁺, and one Cu²⁺.

The aim of studying Ca-MMT was to provide a reference for comparison with Cu-MMT, due to their shared divalent cation nature. Additionally, the difference in ionic sizes between Ca²⁺ and Cu²⁺ provides insights into how cation size affects the swelling pressure and interaction energies in MMT. This comparative analysis was intended to improve the accuracy of MD calculations for Cu-MMT. The inclusion of Ca²⁺ and Na⁺ ions within the original ClayFF package facilitated a more precise evaluation, leveraging the established parametrizations for these ions.
Convergence tests were performed, revealing that a kinetic energy cut-off for wave functions of 1224.5 eV (90 Ry) for Cu-MMT and 680.28 eV (50 Ry) for Ca-MMT, along with a 2 × 2 × 2 Monkhorst-Pack k-point grid [51], were sufficient. Next, geometric optimization of Cu- and Ca-MMT was performed. For this purpose, seven distinct Ca-MMT and twelve Cu-MMT configurations, differing in the Ca/Cu positions on the MMT surface, were constructed. The choice of configurations reflects the differing levels of prior research: since Ca-MMT has been extensively studied and its position within MMT is well understood, fewer configurations were necessary. In contrast, more configurations were used for Cu-MMT to thoroughly explore and find the most stable configuration A periodic simulation box with an initial dimension of 11 Å along the Z axis was employed. Volume and ion relaxation were run for each configuration using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [52]. The total energy and force convergence thresholds were set at 1 × 10⁻⁴ eV and 1 × 10⁻³ atomic units of force (where 1 atomic unit of force is e²/a₀²), respectively. In our DFT setup, we did not explicitly define the charge for the MMT and Ca/Cu. Instead, we used a single Cu or Ca atom. To determine the optimal configuration, the adsorption energy of the adsorbed atoms was computed using Eq. (1):

$$E_{ads} = E_{Ca/Cu-MMT} - (E_{Ca/Cu} - E_{MMT}). \quad (1)$$

Here, $E_{Ca/Cu-MMT}$ represents the energy of the Ca/Cu-MMT system after optimization, $E_{Ca/Cu}$ denotes the energy of a single Ca/Cu atom in a simulation box and $E_{MMT}$ signifies the energy of the optimized MMT structure.

## 2.2. Interatomic interaction potential

This subsection introduces ClayFF, which serves as the foundation for our MD simulations. Then, we detail the process of extending ClayFF for Cu²⁺ interactions.

### 2.2.1. ClayFF

ClayFF [23,38] is constructed using mostly non-bonded interactions including van der Waals in the form of Lennard-Jones (12-6) (L-J) function Eq. (2) and Coulombic interaction energies Eq. (3):

$$E_{VDW} = \sum_{i \neq j} 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right]. \quad (2)$$

$$E_{Cout} = \frac{e^2}{4\pi\epsilon_0} \sum_{i \neq j} \frac{q_i q_j}{r_{ij}}. \quad (3)$$

The parameters $\epsilon_{ij}$ and $\sigma_{ij}$ denote the potential energy well depth and the distance at which the potential energy crosses 0. For simplicity, these parameters are presented for identical atoms, referred to as diagonal interaction terms. Off-diagonal terms are calculated by the Lorentz-Berthelot combination rules which use arithmetic and geometric means for $\sigma_{ij}$ and $\epsilon_{ij}$ [53], respectively. ClayFF includes the flexible simple point charge (SPC) water model introduced by Berendsen et al. [54]. Within this model, each atom is represented by a partial charge at the center of the atom. Specifically, oxygen and hydrogen atoms within water molecules carry partial charges of -0.82 e and +0.41 e, respectively. While ClayFF maintains similar L-J parameters for all oxygen atoms as introduced in the SPC water model, the partial charges vary based on local coordination and environments. ClayFF was originally designed for infinite (periodic) clay sheets; however, several studies have modified ClayFF to investigate clay particles containing hydroxyl edges [23,39,55-60] by extending harmonic or Morse potential for O-H bonds, as well as harmonic potentials for the metal-O-H angles. In our present investigation, a harmonic potential is employed for the bonds and angles.

### 2.2.2. Training procedure for Cu²⁺ interactions

In order to train and validate a ClayFF Cu²⁺ interaction model, molecular statics and MD simulations were conducted in the large-scale atomic/molecular massively parallel simulator (LAMMPS) software package [61,62], employing the modified ClayFF force-field [23]. For all simulations, a cut-off radius of 9.5 Å along with Ewald summation to compute long-ranged electrostatic interactions in reciprocal space was applied.

To initiate the process, we employed a periodic 2 × 1 × 1 unit cell of MMT mimicking the setups utilized in DFT calculations. Within this unit cell, we considered twelve distinct Cu-MMT setups, mirroring those employed in DFT. For each of these unique Cu-MMT configurations, we performed energy minimization via a conjugate gradient algorithm [63] with at least 14 different sets of Cu²⁺ L-J parameters sourced from existing literature [64-77], along with 40 sets of L-J parameters in a range of $\epsilon$ from 0.001 to 0.9 kcal/mol and $\sigma$ from 2.11 to 3.51 Å. This range was selected based on a set of L-J parameters found in the literature that gave us the closest structural similarity to the DFT results. Our aim here was to identify a set of L-J parameters that best aligned with the geometric outcomes obtained from DFT.

Having identified a set of parameters that provided geometry, specifically matching the d-spacing and Cu-O, Cu-Si distances, similar to our DFT results, we proceed with validation of the extended ClayFF.

First, we evaluated the energy of each configuration in two contexts:

1. The energy difference of each configuration in ClayFF relative to the most stable structure identified in ClayFF.

2. The energy difference of each configuration in DFT relative to the most stable structure identified in DFT.

We considered 12 configurations for Cu-MMT and 7 configurations for Ca-MMT by varying the positions of ions, and then calculated the energy differences to determine the stability of each configuration. To achieve this, we optimized each structure with both DFT and the extended ClayFF model. This allowed us to identify the local minima for each initial configuration and subsequently compare the energy of each configuration in relation to the lowest found energy in either DFT or ClayFF. This approach was implemented for both Ca- and Cu-MMTs to ensure a comprehensive assessment of the parameterization.

Second, we assessed the spatial arrangement of individual atoms within the Cu-MMT system as determined by ClayFF and compared it to that from DFT calculations. This validation procedure was applied to the Ca-MMT system. Our principal objective in this validation was to ascertain that the interatomic interactions involving Cu²⁺ within the ClayFF force-field had achieved a level of consistency on par with that attained for Ca²⁺ within the original ClayFF force-field.

In another validation test, since bentonite is a swelling clay, we conducted separate simulations with similar MMT structures containing three water molecules, one Ca²⁺, one Cu²⁺, and two Na⁺ ions in the canonical (NVT) ensemble at 500 K using the Nose-Hoover thermostat. Forty configurations for each structure were extracted from the dump file and examined in SCF DFT calculations to compare the correlation between atomic forces generated by DFT, the original ClayFF for Na⁺ and Ca²⁺, and the extended ClayFF for Cu²⁺.

Finally to investigate the behavior of aqueous Cu²⁺, we initially created a large simulation box containing 2790 SPC water molecules with periodic boundary conditions in all dimensions. We simulated it in an isothermal-isobaric (NPT) ensemble at 300 K and 1 atm pressure. We assessed equilibration by monitoring density and the radial distribution function (RDF) of H-H, H-O, and O-O. The average density of equilibrated water during a 500 ps period was found to be 0.994 g/cm³, consistent with values from other sources [78]. Afterward, one Cu²⁺ was inserted in the middle of the simulation box, along with two Cl⁻ placed at the farthest distance from Cu²⁺ to minimize their interactions. The reason behind having two Cl⁻ was that the system had to be neutral in order to employ the Ewald summation. Water molecules within a 2 Å radius of Cu²⁺ were removed. Subsequently,
the system underwent three separate simulations, each focusing on a single set of L-J parameter. Each system was pre-equilibrated in an NPT ensemble at 300 K and 1 atm for 100 ps, then equilibrated quantities were sampled in the canonical (NVT) ensemble at 300 K using Nose-Hoover thermostat for a duration of 500 ps. Validation of the new force fields parameters involved assessing geometric characteristics of the Cu²⁺ ion and its surroundings in water, including RDF, coordination number (CN) and angular distribution function (ADF).

## 2.3. Potentials of mean force

To gain insights into the behavior of MMT in the presence of Na⁺, Ca²⁺, and Cu²⁺ ions at the atomic scale, potentials of mean force (PMF) were computed. A PMF quantifies the variation in free energy between different states along a reaction coordinate, connecting the energies of the states of interest through a defined pathway [79–81]. In our work, PMFs are derived using force constraint technique [81], focusing on edge-to-edge and center-to-center MMT configurations. This technique involves sampling the forces acting on fixed reaction coordinates, as described by Eq. (4):

$$ \Delta A_{a-b} = \int_{r_a}^{r_b} F_r dr. \qquad (4) $$

Here, $\Delta A_{a-b}$ represents the difference in free energy between states ‘a’ and ‘b’, and it is computed by integrating the force required to maintain the MMT at a specified distance of $r_a$ and $r_b$.

Three different models of MMTs were built using the ATOM package [50] to calculate PMFs.

1. A small hexagonal MMT platelet which aligns with the methodology detailed in Zhang et al. work [27]. The prior research by Underwood and Bourg emphasized the energetic favorability of hexagonal-shaped MMTs due to their minimized number of broken silicate rings [55,59]. This model maintains a stoichiometric ratio of 4.25:1 for Al:Mg with the charge of -8. The structural formula for the Cu-MMT is

[Cu₄].[Si₁₆₈·[Mg₁₆Al₆₈]O₁₉₂(OH)₇₈].

2. A periodic sheet of MMT with the structural formula

[Cu₁₃·[Si₃₂₀·[Mg₂₆Al₁₃₄]O₈₀₀(OH)₁₆₀].

3. A semi-periodic sheet, non-periodic in the x and periodic in y direction, with the structural formula

[Cu₁₆·[Si₃₃₆·[Mg₂₈Al₁₄₀]O₈₂₄(OH)₂₀₀].

In the small platelet and semi-periodic model hydroxyl groups were introduced to rectify the edges. Al-O-H and Si-O-H in the edges were integrated into the simulation via the extension of harmonic bonds and angles from the initial ClayFF force field [58]. Initially, a large box of water was structurally optimized using the conjugate gradient (CG) method. Subsequently, equilibration occurred under an NPT ensemble at a temperature of 300 K and a pressure of 1 atm to ensure proper thermal, pressure, and energy stabilization of the system. Following this, Na-MMT, Ca-MMT, and Cu-MMT were each introduced into the box of water through separate simulations. Each PMF graph was generated through a sequence of simulations, wherein the distances between two MMT platelets or sheets were incrementally increased by an average of 0.2 Å. In each scenario, water molecules within a 2 Å radius of each platelet were excluded to prevent overlapping water molecules with MMT platelets and to avoid undesired interactions in the initial steps of simulation.

Across all simulations, metals in the octahedral layers remained fixed in all x, y, and z directions: 11 out of 42 for platelets, 12 out of 160 for periodic sheets, and 14 out of 168 for semi-periodic sheets.

Subsequently, every system underwent pre-equilibration under NPT conditions at 300 K and 1 atm for 50 ps, followed by further equilibration and force sampling under NVT conditions at 300 K for 1 ns. A timestep of 1 fs was chosen for all simulations.

<caption>**Table 1**</caption>

Structural details of optimized Cu-MMT in DFT and our best-fitted L-J parameters ($\sigma = 2.13$ Å) and $\epsilon = 0.05$ kcal/mol) for Cu²⁺ in ClayFF.

<table><thead><tr><th>Method</th><th>d-spacing (Å)</th><th>Cu-Si distance (Å)</th><th>Cu-O distance (Å)</th></tr></thead><tbody><tr><td>DFT</td><td>9.54</td><td>2.75</td><td>2.08, 2.06</td></tr><tr><td>ClayFF ($\sigma = 2.13$ and $\epsilon = 0.05$)</td><td>9.35</td><td>2.67</td><td>2.02, 2.02</td></tr><tr><td>ClayFF ($\sigma = 2.16$ and $\epsilon = 0.05$)</td><td>9.50</td><td>2.66</td><td>2.05, 2.02</td></tr><tr><td>ClayFF ($\sigma = 2.21$ and $\epsilon = 0.05$)</td><td>9.50</td><td>2.68</td><td>2.06, 2.05</td></tr></tbody></table>

Swelling pressure calculations for Na-MMT, Ca-MMT, and Cu-MMT were determined by analyzing the gradients of PMFs at various dry densities. These dry densities, defined as $\frac{m}{A \times d}$, where m is the mass of the clay with ions, A is the area of the clay, and d is the d-spacing, were calculated based on the center-to-center distances between two cation-MMTs, as inspired by the approach detailed in Sun's study [41]. In addition to this analysis, we also computed the mean square displacement (MSD) for each counter-ion within the MMT clays.

# 3. Results and analysis

## 3.1. Cu²⁺ interatomic interactions model

### 3.1.1. Results of the fitting procedure

We start by presenting our DFT-calculated adsorption energies, providing an insight into the strength of interaction between cations and MMT surface. An adsorption energy of -255.73 kcal/mol for Cu²⁺ cations was obtained. The optimal position for Cu²⁺ was found just above the MMT surface, nearly in the middle of interlayer region, positioned close to the boundary of the ring, directly above one of the two Al→Mg substitutions in the octahedral layer. Cu²⁺ exhibited a bonding distance of 2.75 Å with Si atom and 2.08 Å and 2.06 Å with O atoms surrounding the Si atom in the tetrahedral layer of MMT. Details provided in Table 1, and Fig. 1 illustrates the corresponding Cu-MMT structure. An important observation here is the distinct behavior of Cu²⁺ compared to Ca²⁺ in terms of adsorption on MMT surfaces. The adsorption energy for Ca²⁺ for the corresponding configuration shown in Fig. 1 was determined to be -147.82 kcal/mol. Interestingly, in the case of Ca²⁺, the ion is positioned precisely in the center of a ring and in close proximity to one of the MMT surfaces, consistent with prior calculations [82]. By comparing the structure of Cu²⁺ shown in Fig. 1 with the calculated structure for Ca²⁺, where the Ca²⁺ ion is positioned in the center of the tetrahedral ring of MMT, we can understand the differences in adsorption characteristics between Cu²⁺ and Ca²⁺.

After identifying the most stable Cu-MMT configuration using DFT, we aimed to replicate this structure with ClayFF. We performed a search involving 14 parameter sets to achieve accurate d-spacing, Cu-Si and Cu-O distances. The best-fitting L-J parameters were determined to be $\sigma = 2.13$ Å and $\epsilon = 0.05$ kcal/mol, as were proposed by K. Anitha et al. [66]. To further enhance the structural accuracy, we explored variations in $\sigma$, including $\sigma = 2.16$ and 2.21 Å, as detailed in Table 1. $\sigma = 2.21$ Å resulted in a structure closely resembling the DFT results. Fig. 1(b) illustrates the optimized Cu-MMT structure associated with the extended ClayFF. Note that Cu²⁺ as described by ClayFF resides in a different ring as compared to DFT, while maintaining similar bonding lengths and d-spacing as observed in the DFT simulations. In an effort to validate the positioning of Cu²⁺ ions, we attempted to place the Cu²⁺ ion near the DFT-derived position and subsequently relaxed the structure. However, during the relaxation process, the Cu²⁺ ion consistently migrated to another ring and did not remain at the initial DFT position.

In addition to the parameter sets explored in this study, we evaluated a set of parameters reported by Zhang et al. in a recent paper [76]:
<figure><img src="image_4.png" /><figcaption>Fig. 1. Structural comparison between (a) DFT-optimized periodic Cu-MMT and (b) Cu-MMT optimized using the extended ClayFF with $σ = 2.21$ Å and $ε = 0.05$ kcal/mol. The arrows indicate the d-spacing values. In this representation, brown, red, green, gray and beige spheres represent Cu$^{2+}$ ions, oxygen, magnesium, aluminum, and silicon atoms, respectively. The right side figures show a top view of the structure with the bottom tetrahedral and octahedral layers made transparent for better illustration.</figcaption></figure>

$σ = 3.114$ Å and $ε = 0.005$ kcal/mol. This L-J model closely approximated our findings, with Cu-Si distances at 2.67 Å and Cu-O distances at 2.14 Å and 2.09 Å, along with a d-spacing of 9.64 Å. While these parameters effectively described the Cu-MMT system in a dry environment, similar to our DFT results, our developed model exhibited higher accuracy. It is worth noting that the value for $ε$ is expressed in a different unit compared to the cited value in Ref. [76], but they are equivalent.

We conducted a series of calculations to assess how Cu-MMT energies vary as Cu$^{2+}$ traverses the interlayer space. This involved testing 40 sets of L-J parameters ranging from $σ = 2.11$ to 2.91 Å and $ε = 0.001$ to 1 kcal/mol, in addition to the parameters employed by Zhang et al. [76]. The results, illustrated in Fig. 2, demonstrate the outcomes obtained using Zhang's parameter set.

The results, illustrated in Fig. 2, demonstrate the outcomes obtained using Zhang's parameter set compared with our results and the DFT results. The zoomed-in graph highlights the changes in energy around the minimum well, where the probability of finding Cu$^{2+}$ is highest. Even under extreme stresses, the Cu$^{2+}$ ion remains within the range of 9 to 9.5 Å from its equilibrium position. Energy curves beyond this range do not significantly impact the mechanical properties due to the large energy scale compared to thermal agitation energy (kT). Our model is similar to Zhang's, both depicting similar potential energy changes at the thermal scale with comparable well depths. The blue line matches the DFT results around the minimum as accurately as Zhang's model, indicating similar interaction characteristics.

### 3.1.2. Validation

First, the positions of each element in the ClayFF-optimized Cu-MMT and their counterparts in DFT simulations were compared. The

atomic coordinates from both simulations were converted to fractional space for accurate comparison. One atom was used as a reference in both the DFT and ClayFF simulations, and the positions of all other atoms were normalized relative to this reference atom. The deviations in the position of each element between ClayFF and DFT, normalized to the unit cell dimensions, are visually presented in Fig. 3. We also optimized Ca-MMT using both DFT and ClayFF and quantified the deviations in atomic positions from the DFT-optimized structure. The extended Cu$^{2+}$ ClayFF model led to deviations for the MMT atoms comparable to those associated with the Ca$^{2+}$ interaction model in ClayFF inside MMT. Although the deviation for Cu$^{2+}$ is higher than for Ca$^{2+}$, the bonding length and d-spacing are similar to DFT. This higher deviation for Cu$^{2+}$ can be attributed to the movement of Cu to another ring, which caused relatively high deviation for a small cell size.

Another analysis involved comparing the Ca$^{2+}$ and Cu$^{2+}$ adsorption energies for each configuration relative to the minimum energy state using either ClayFF or DFT. For this analysis, we optimized each configuration to find the local minimum around the initial configuration and compared the obtained energies. In Fig. 4, each point corresponds to an initial configuration. As illustrated, the extended ClayFF successfully described the lowest-energy Cu$^{2+}$ site. Notably, while similar initial configurations led to the lowest energy for both DFT and ClayFF, in ClayFF, Cu$^{2+}$ moves to a different ring than in DFT. Specifically, in DFT optimization, Cu$^{2+}$ positions itself above the Al→Mg substitutions in the octahedral layer, whereas in ClayFF, the ion moves above an Al atom in a ring without substitution. However, the distances between Cu$^{2+}$ and O and Si atoms in the tetrahedral layer of MMT are very similar in both methods. The DFT-ClayFF correlation pertaining to Cu-MMT is more robust than that pertaining to Ca-MMT. According to DFT, the variation of the adsorption energies in Cu-MMT are larger than those
<figure><img src="image_5.png" /><figcaption>Fig. 2. A comparative analysis of the Cu-MMT energy profile along the perpendicular direction of MMT surfaces. The energy profiles are generated from DFT (black line), our best-fitted L-J parameters (blue line), and the parameters proposed by Zhang et al. (red line) [76].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Comparison of particle positions as predicted by ClayFF and DFT after structural relaxation. The distances are normalized to the unit cell dimensions. Ca-MMT and Cu-MMT deviations are shown in blue and red, respectively, with overlapping deviations displayed in purple.</figcaption></figure>

involved in Ca-MMT. ClayFF captures this effect, albeit exaggeratedly so. DFT suggests Cu²⁺ adsorption energies range between 0 and 17 kcal/mol-relative to the lowest-energy site–ClayFF suggest a 0 to 35 kcal/mol range, as compared to 0 to 8 kcal/mol for Ca-MMT. Notably, low-energy Cu-MMT configurations, which are of higher importance, show relatively good correlation between DFT and ClayFF. However, low-energy configurations of Ca-MMT in DFT (energies around 0.1 kcal/mol) are not well predicted by ClayFF.

To validate the force field, we conducted NVT ensemble simulations at 500 K for Na-, Ca-, and Cu-MMT structures with three water molecules using ClayFF. We extracted forty configurations for each structure and performed SCF DFT calculations to compare the correlation between atomic forces generated by DFT and ClayFF (original for Na⁺ and Ca²⁺, and extended for Cu²⁺). The resulting forces on atoms inside MMT, including water molecules, are shown in Fig. 5, and the forces on the counter-ions are shown in Fig. 6. Each point represents a force calculated by SCF DFT (x-axis) and ClayFF (y-axis). The solid line indicates the linear regression fit, while the dashed line represents the ideal 1:1 correlation. Fig. 5 shows the comparison of atomic forces on atoms inside the Na- (green), Ca- (blue), and Cu-MMT (red) systems.

For the Na-MMT system, the data points are closely clustered around the 1:1 line, indicating strong agreement between DFT and ClayFF forces. The high R² values and low RMSE demonstrate the accuracy of the original ClayFF for Na⁺ ions. In the Ca-MMT system, which is included in the original ClayFF, the results show good agreement with DFT, albeit with higher RMSE and lower R² values compared to Na-MMT. This reflects the inherent reliability of ClayFF for Ca²⁺ ions. Interestingly, the Cu-MMT system exhibits statistical trends similar to the Ca-MMT system. The slopes and intercepts for Cu-MMT are close to those of Ca-MMT, with comparable R² values and RMSE. This suggests that the extended ClayFF for Cu²⁺ performs reasonably well, achieving a level of accuracy similar to Ca²⁺, which is originally parameterized in ClayFF.

In Fig. 6, we observe the comparison of atomic forces on ions in the Na- (green), Ca- (blue), and Cu-MMT (red) systems. For the Na-MMT system, the forces on Na⁺ ions show a high degree of correlation with the DFT results, evidenced by slopes and intercepts close to the ideal values, high R² values, and low RMSE. This indicates a strong agreement between ClayFF and DFT forces. In the Ca-MMT system, the comparison shows good agreement with DFT forces for Ca²⁺ ions. The
<figure><img src="image_7.png" /><figcaption>Fig. 4. Adsorption energy – setting the lowest-energy site at zero – comparison. ClayFF and DFT calculations describing (a) Cu²⁺ and (b) Ca²⁺ adsorbed on different MMT surface sites after structural relaxation are compared. Dashed lines represent a one-to-one relationship.</figcaption></figure>

slopes and intercepts, although slightly different from the ideal values, still indicate a high correlation, with reasonably high R² values and low RMSE, though not as strong as in the Na⁺ system. For the Cu²⁺ ions, while there is more deviation compared to Na⁺, the statistical measures (slopes, intercepts, R² values, and RMSE) in the y-direction are similar to those of the Ca²⁺ ions, and results in other directions show similar accuracy to those for atoms inside Ca-MMT. It should be noted that Cu²⁺'s partially filled d-orbitals lead to more complex chemical properties and structures compared to Ca²⁺ and Na⁺. Consequently, performing classical force field calculations for Cu²⁺, such as the simple L-J interaction energies used in ClayFF, does not achieve the same level of accuracy as for simpler ions like Ca²⁺ and Na⁺.

Validation of the extended ClayFF model in a hydrated environment involved assessing the geometric properties of Cu²⁺ ions immersed in water, including RDF, CN and ADF. We considered the three best-performing sets of L-J parameters found as detailed in Table 2. Fig. 7 highlights the results achieved using the optimal parameters, σ = 2.21 Å and ε = 0.05 kcal/mol. In Fig. 7, the RDF for Cu-O<sub>water</sub> is depicted by the blue line, while the hydration number of Cu²⁺ is shown in red. The first peak in the RDF corresponds to the average distance between Cu²⁺ and the oxygen atoms in its first hydration shell, measuring 2.04 Å. A visual representation of Cu²⁺'s first hydration shell is provided in Fig. 7. The second peak in the RDF, located at 4.24 Å, signifies the distance between Cu²⁺ and the oxygens in the second hydration shell. Coordination numbers were determined by integrating the Cu-O<sub>water</sub> RDF. The analysis of the first coordination number, representing the count of oxygen atoms in the first hydration shell, confirms that the new force field favors a 6-fold coordination for Cu²⁺ cations in the water solution, consistent with previous studies, as elaborated in Table 2 [83–87]. The second-shell coordination number, calculated for distances between 3 Å and 4.8 Å, was found to be 12.3, closely aligning with results from quantum mechanical/molecular mechanical molecular dynamics (QM/MM MD) simulations [85], which are also compared in Table 2. Additionally, the ADF of O-Cu-O for the first hydration shell around the Cu²⁺ ion, as presented in Fig. 7, indicated angles of 89.5° and 175.5°. These results demonstrate excellent agreement with findings from other studies [85–90]. However, note that certain quantum effects observed using QM approaches, such as Jahn-Teller inversion, cannot be replicated using classical force-fields.

## 3.2. Potentials of mean force

We conducted a series of MD simulations using the extended ClayFF force field to calculate potentials of mean force (PMF) profiles in wet MMT systems. Our study includes four different scenarios, including both center-to-center and edge-to-edge configurations of MMT

platelets. The rationale behind exploring various systems lies in the recognition that MMT platelets in real-world settings exhibit a wide range of sizes, shapes, and spatial arrangements when interacting with each other.

### 3.2.1. Measured potentials of mean force

Our investigation was initiated by examining the behavior of two center-to-center MMT platelets, including Na⁺, Ca²⁺, and Cu²⁺, in separate simulations. The resulting average PMF curves, along with a snapshot illustrating the Na-MMT structure at a specific center-to-center separation, are presented in Fig. 8. Remarkably, the PMF curves, irrespective of the cation type, reveal a rich landscape of local energy minima, corresponding to various hydration states. These states encompass configurations with zero water layers (0-W), one water layer (1-W), and two water layers (2-W) inserted between the MMT platelets. Notably, the transitions between the 0-W and 1-W states, as well as between the 1-W and 2-W states, consistently occur at distances of approximately 2.5–3 Å. This distance closely matches the diameter of a water molecule and aligns with experimental measurements of hydration forces observed between diverse surfaces, including mica and MMT materials [34,35,93,94]. There is a pronounced repulsive interaction at distances less than 16 Å. There are substantial energy barriers situated at 10.5–11.5 Å for the hydration process from the 0-W to 1-W state. These energy barriers can be partly attributed to the contribution of hydrogen bonding at the edges of the platelets. This observation is consistent with a prior study [39], which demonstrated that hydrogen bonds at the platelet edges serve as gatekeepers, controlling the entry of water molecules into the interlayer spaces of MMT. This effect primarily influences the first layer of hydration. The initiation of the swelling process is found to depend on breaking hydrogen bonds and was influenced by the size and shape of the MMT platelets [39]. The differences between the two local energy minima indicate the relative stability of different hydrated states. The significant difference between minima at 0-W and 1-W for all structures underscores the strong preference of MMT platelets for water absorption. This behavior is expected for MMT in contact with liquid water under standard conditions, as opposed to low hydration states which become more favorable in vapor when relative humidity (RH) is reduced. Furthermore, the energy difference between 0-W and 1-W of Na-MMT is lower than that observed for Ca- and Cu-MMT, possibly due to variations in the hydration energy of the first water shell around different cations. Additionally, the larger size of hydrated Cu²⁺ leads to its local energy minima being located at a slightly greater distance as compared to the other two cations.

Additionally, to address concerns about the anisotropy of platelets and their rotational degrees of freedom, we included PMF calculations for center-to-center interactions with a 30° rotation applied to one of
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Comparison of atomic forces on atoms inside the Na- (green), Ca-(blue), and Cu-MMT (red) systems, calculated using SCF DFT (x-axis) and ClayFF (y-axis). Each row corresponds to a different ion type, while the columns represent forces in the x, y, and z directions respectively. The slope, intercept, their 95% confidence interval, R² value, and RMSE are provided for each fit. The solid line indicates the linear regression fit, and the dashed line represents the ideal 1:1 correlation.</figcaption></figure>

<caption>**Table 2**</caption>

Comparison of Cu²⁺ hydration shell properties obtained using various methodologies, including Cu-O<sub>Water</sub> distances and coordination numbers (CN) for the first and second hydration shells, along with the O-Cu-O angle in the first hydration shell.

<table><thead><tr><th>Method</th><th>Cu-O<sub>1st</sub> shell (Å)</th><th>CN<sub>1st</sub></th><th>Cu-O<sub>2nd</sub> shell (Å)</th><th>CN<sub>2nd</sub></th><th>Angle<sub>O-Cu-O</sub> (°)</th></tr></thead><tbody><tr><td>ClayFF (σ = 2.13, ε = 0.05)</td><td>2.02</td><td>6</td><td>4.28</td><td>12.8</td><td>88.5, 175.5</td></tr><tr><td>ClayFF (σ = 2.16, ε = 0.05)</td><td>2.02</td><td>6</td><td>4.27</td><td>12.64</td><td>88.5, 174.5</td></tr><tr><td>ClayFF (σ = 2.21, ε = 0.05)</td><td>2.04</td><td>6</td><td>4.24</td><td>12.3</td><td>89.5, 175.5</td></tr><tr><td>HF/MM+3bd MD [83]</td><td>2.08</td><td>6</td><td>4.2</td><td>14.5</td><td></td></tr><tr><td>EXAFS/LAXS [84]</td><td>1.95, 2.29</td><td>6</td><td>4.17</td><td>8</td><td></td></tr><tr><td>B3LYP QM/MM [85]</td><td>2.02, 2.29</td><td>6</td><td>4.13</td><td>11.9</td><td>89, 173</td></tr><tr><td>CPMD [91]</td><td>1.96, 1.96</td><td>5</td><td></td><td></td><td></td></tr><tr><td>CPMD [88]</td><td>2, 2.45</td><td>5</td><td>4.03</td><td>8</td><td>90, 180</td></tr><tr><td>ReaxFF [86]</td><td>1.94, 2.27</td><td>6</td><td>4.27</td><td>12.5</td><td>90, 176</td></tr><tr><td>Neutron diffraction [87,92]</td><td>1.97</td><td>6</td><td></td><td></td><td></td></tr><tr><td>X-ray diffraction [87,89]</td><td>1.98, 2.39</td><td>6</td><td>3.96</td><td>11.1</td><td></td></tr><tr><td>X-ray diffraction [87,89]</td><td>1.98, 2.34</td><td>6</td><td>3.95</td><td>11.6</td><td></td></tr><tr><td>X-ray diffraction [87,90]</td><td>2.01, 2.33</td><td>6</td><td>4.2</td><td>7.6</td><td></td></tr></tbody></table>

the platelet surfaces, as shown in Fig. 8(b). The general trend remains the same, and hydrated layers appear at similar d-spacing. However, the barrier energy from 0-W to 1-W has increased in Ca- and Cu-MMT, possibly due to the more pronounced effect of hydrogen bonding at the edges. The difference in PMF profiles before and after rotation, especially for the Cu-MMT structure, may be attributed to the specific interactions of Cu²⁺ ions with the bridge oxygen atoms, leading to distinct binding behaviors that are sensitive to the orientation of the platelets. It is worth mentioning that while a 0-W state is included in our simulations for theoretical completeness, in practical scenarios, some level of hydration is always present. However, in vapor at very low relative humidity (RH) or high temperatures, low hydration states can indeed become more favorable.

In our study, we conducted PMF calculations for interactions involving two edge-to-edge platelets according to Fig. 8(c). This type

of interaction does not exhibit significant dependence on the cations present in the environment. While the local minima in edge-to-edge interactions are not as pronounced as those observed in center-to-center interactions, we observe an oscillation pattern in the PMF profiles. This pattern is attributed to the presence of water layers positioned between the two platelets.

Another assessment involved studying the interaction between a single small platelet and an infinite (periodic) sheet to replicate the interaction between small and significantly larger platelets. The resulting PMF and an illustration of the system are presented in Fig. 9. The average PMF curves share similarities with those of individual platelets. However, the initial energy barrier for the transition from 0-W to 1-W is lower compared to that observed for small platelets. This discrepancy can be attributed to the reduced contribution of hydrogen bonds in these interactions. The interaction energies exhibit oscillatory patterns
<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison of atomic forces on ions in the Na- (green), Ca- (blue), and Cu-MMT (red) systems, calculated using SCF DFT (x-axis) and ClayFF (y-axis). Each row corresponds to a different ion type, while the columns represent forces in the x, y, and z directions, respectively. The slope, intercept, their 95% confidence interval, R² value, and RMSE are provided for each fit. The solid line indicates the linear regression fit, and the dashed line represents the ideal 1:1 correlation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) Radial distribution function (RDF) of Cu-O<sub>water</sub> for hydrated Cu<sup>2+</sup> in the flexible SPC water model using the extended ClayFF force field. The solid blue line represents the RDF and the dashed red line represents the coordination numbers (CN). (b) Angular distribution function (ADF) of O<sub>water</sub>-Cu-O<sub>water</sub> in the first hydration shell using the extended ClayFF force field.</figcaption></figure>

as a consequence of the hydration forces that come into play when
increasing the size of a single MMT layer, while also decreasing the
influence of hydrogen bonds.

The latest set of calculated PMF profiles, as depicted in Fig. 10,
pertains to the interaction between two semi-infinite (semi-periodic)
sheets of MMT. This system effectively mimics the behavior of very
large platelets while considering the influence of hydrogen bonds at the
edges of these sheets, providing valuable insights into the interaction
dynamics at this larger scale. As illustrated in Fig. 10, the overall
shape of the PMF curves exhibits similarities with those observed for

smaller platelets. One important feature, as revealed by the zoomed-
in PMF graphs in Fig. 10, is the distinct behavior regarding water
adsorption when the d-spacing exceeds approximately ~ 21 Å. In this
specific region, Na-MMT continuously adsorbs water without encoun-
tering any significant energy barriers. It appears to transition beyond
the crystalline region, a phenomenon that aligns with a previous study
on Arizona Na-MMT, which demonstrated the absence of forbidden
layers in its structure [26]. In contrast, Ca-MMT and Cu-MMT face
relatively minor energy barriers in this d-spacing range and do not
exhibit swelling. This disparity in behavior could be attributed to the
presence of forbidden layers, as observed in the work of Meleshyn et al.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Potentials of mean force (PMF) for (a) center-to-center interactions of Na-, Ca-, and Cu-MMT platelets, (b) center-to-center interactions with a 30° rotation applied to one of the platelet surfaces and (c) edge-to-edge interactions of the aforementioned structures. The PMF results, represented by dots, are averaged, and the shaded colors indicate the variability obtained from 100 bootstrap samples.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Potentials of mean force (PMF) for center-to-center configuration between one periodic sheet and one small platelet of Na-, Ca- and Cu-MMT. (a) The PMFs, represented by dots, are averaged, and the shaded colors indicate the variability obtained from 100 bootstrap samples. The zoomed-in graph reveals a well-defined third water minimum. (b) Corresponding structure for Na-MMT is shown.</figcaption></figure>

for Wyoming Na-MMT [32]. The concept of forbidden layers arises from the chainlike structure formed by Na⁺ ions and water molecules in the interlayer region. In the case of Wyoming Na-MMT, Na⁺ ions form an inner complex, however; in Arizona MMT substitution occurs in the octahedral layer and Na⁺ forms outer-sphere complex. This leads to the absence of forbidden layers in Arizona Na-MMT. However, it is plausible that Ca²⁺ and Cu²⁺ ions, due to their higher charges, form inner complexes with the Arizona MMT layers, resulting in the creation of forbidden layers. An experimental study using low-angle XRD has demonstrated that Ca-Wyoming, Ca-Otay bentonites, and Ca-hectorite undergo a stepwise expansion, ultimately expansion stops at d-spacing of 19.6 Å for Ca-Wyoming and 22 Å for other two clays [26]. These findings support the presence of forbidden layers in Ca-MMT.

### 3.2.2. Swelling pressure

In our study, we aimed to quantify and compare the swelling pressures induced by Cu-MMT systems with those of Na-MMT and Ca-MMT. To the best of our knowledge, limited information is available regarding the swelling pressures of Cu-MMT.

We employed the slopes derived from average PMFs to interpolate swelling pressures of Na-, Ca-, and Cu-MMT. The derivatives were obtained using cubic spline interpolation and the Savitzky-Golay filter. Additionally, we determined distances to establish dry densities, relying on the methodology introduced by Linlin Sun et al. [41],

which involved the conversion of d-spacing data into dry densities. As previously defined, these dry densities are calculated using Eq. (5):

$$ \frac{m}{A \times d'} \qquad (5) $$

where m is the mass of the clay with ions, A is the area of the clay, and d' is the d-spacing. The PMF is equivalent to the grand potential in our system and can be expressed as Eq. (6):

$$ \mathrm{PMF}(h) = \int_{h_0}^{h} P(h') \, dh'. \qquad (6) $$

At controlled pressure (P<sub>control</sub>), the minimized potential is the osmotic potential (λ):

$$ \lambda = \mathrm{Grand\ Potential} + P_{\mathrm{control}} \cdot h. \qquad (7) $$

In thermodynamic equilibrium, this osmotic potential is minimized. The difference in osmotic potential between two states can be expressed as Eq. (8):

$$ \lambda - \lambda_0 = \int_{h_0}^{h} (P - P_{\mathrm{control}}) \, dh. \qquad (8) $$

The value of P<sub>control</sub> that provides the transition point is found where this integral equals zero. This indicates that the areas under and above the controlled pressure line on the swelling pressure curve are equal, ensuring phase equilibrium and indicating a phase transition [95,96].

The obtained swelling pressures from our simulations are visually presented in Fig. 11, marked with filled circle symbols and lines. Solid lines indicate stable regions, dashed lines indicate metastable regions
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Potentials of mean force (PMF) for center-to-center orientation of two semi-periodic (a) Na-, (b) Ca- and (c) Cu-MMT. The PMF results, represented by dots, are averaged, and the shaded colors indicate the variability obtained from 100 bootstrap samples. Interestingly, the zoomed-in graph shows a higher slope for the Na-MMT interaction compared to Ca and Cu-MMT.</figcaption></figure>

and transparent lines indicate unstable regions. The points show the transitions, which are obtained using the explained equations.

The overall trend suggests that the stepwise stable regions can be described by two confining exponential equations. This is consistent with patterns observed in previous experimental and simulation studies with pure water and saline solutions [41,42,97,98].

To validate our simulation results, we compared them with data from experimental studies. In Fig. 11, we include swelling pressure for 24 different bentonites, each with different compositions or exchangeable cations (Clay materials (Exp.1)) [99]. We also reference the swelling pressure data for Wyoming MX-80 bentonite, dominated by either Ca (Ca-MMT (Exp.2)) or Na (Na-MMT (Exp.2)) [100]. Additionally, we present measurements of swelling pressure in bentonite (Exp.3) containing 50%-60% MMT with Ca²⁺ and Mg²⁺ exchangeable cations [98], and compared them with Bucher et al. and Sitz et al. work (Bentonite (Exp. 4-5)) [101,102].

For the entire range of densities, the stability trend of the swelling pressure for all structures lies between two exponential curves. For Na-MMT, the lower exponential curve aligns well with the experimental results, indicated by the empty circled markers [98–102]. This alignment suggests that our model accurately captures the swelling behavior of Na-MMT across this density range. Similarly, for Cu-MMT, the stability trend follows the same pattern. For Ca-MMT, the lower bound of stability aligns with the experimental results for both low and high densities, showing good agreement with experimental data.

At lower densities, below 1.3 g/cm³, Na-MMT exhibits relatively high swelling pressures compared to Ca-MMT and Cu-MMT. As the density increases from 1.3 to 1.65 g/cm³, the swelling pressure of Ca-MMT begins to surpass that of Na-MMT and Cu-MMT. Cu-MMT generally displays intermediate swelling pressures between Na-MMT and Ca-MMT across this density range.

It is important to note that while the lower bound of Na-MMT's swelling pressure aligns with experimental data, the trend observed for Ca-MMT, particularly in the density range of 1.35 to 1.65 g/cm³, where it surpasses the swelling pressure of Na-MMT, contrasts with experimental observations that typically show Na-MMT swelling more than Ca-MMT at these densities. This pattern, however, has been observed in other molecular dynamics simulation studies as well [41, 42]. The deviation from experimental results could be due to several factors, including the limitations of the force field used in the simulations, the assumptions made regarding the ionic environment, and the simplifications inherent in the model.

The overestimation in our calculations could be attributed to the idealized nature of our simulated models, which consist of two perfectly stacked MMT clays. In practice, Effective Montmorillonite Dry Density (EMDD) is used to account for factors influencing swelling mechanisms, and it is typically lower than the actual dry density

(DD). Our comparison with EMDD-adjusted experimental data suggests that while our results are higher, this discrepancy can be attributed to differences in model assumptions and environmental conditions in simulations versus experiments. For instance, the percentage of water used in experiments versus simulations can influence the results, and experimental conditions are highly sensitive to environmental variables such as temperature, humidity, and sample preparation methods. Additionally, the charges on bentonite clay particles can affect the interaction with cations and contribute to discrepancies between simulated and experimental results. Further investigation and refinement of the simulation parameters may be needed to achieve better agreement with experimental data.

### 3.3. Mobility of interlayer cations

The PMF calculations for the semi-periodic MMTs identified approximate local minima for each structure, allowing us to track the trajectories of interlayer ions and subsequently compute their mean square displacements (MSD). After allowing sufficient time for the system to equilibrate, ensuring the amount of water between semi-infinite sheets and the forces reached equilibration, we replaced the semi-infinite platelets with the infinite platelets and measured the MSD based on the average MSD for ions in two interlayer regions.

The MSD graphs presented in Fig. 12 depict the mobility of ions for three different water content scenarios: 1-W, 2-W, and 3-W layers. These graphs show averaged MSD calculations derived from 50 MSD runs over a 5 ns simulation period.

Across all scenarios, Na⁺ ions exhibited higher mobility than the other ions. We quantified ion self-diffusion coefficients (D) and reported them with uncertainties in the form of standard errors, as detailed in Table 3. For the diffusion calculations, we excluded the last 0.1 ns interval to avoid inaccuracies.

The obtained diffusion coefficients revealed that, although an increase in the number of hydrated layers led to a slight initial increase in diffusion in the Z-direction (perpendicular to the surface of MMTs) for all ions, the curves ultimately saturated. This saturation indicates that there was no significant diffusion in this direction across all cases. Therefore, we considered the total diffusion in the XY plane.

The self-diffusion coefficients of Na⁺ ions in the 1-, 2-, and 3-W scenarios are consistent with prior literature, aligning closely with values of 0.53–1.09, 3.6–6.5 and 5.3–7.9 × 10⁻¹⁰ m²/s reported at temperatures of 298 K and 323 K [103].

Similarly, the self-diffusion for Ca²⁺ ions closely match findings from previous MD simulations involving highly charged MMT structures of ~0.75 e per unit cell [104]. This study reported self-diffusion coefficients of 0.6, 1.2, and 2.2 × 10⁻¹⁰ m²/s for 1-, 2-, and 3-W states, respectively, which is relatively consistent with our observations for
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Comparison of simulated swelling pressures of Na+, Ca+, and Cu-MMT with experimental results as a function of dry density [98–102] Solid lines represent stable regions, dashed lines represent metastable regions, and transparent lines represent unstable regions. Filled markers and lines indicate our simulation results, while empty circles represent experimental data. Positive swelling pressures indicate layer expansion due to repulsive forces, while negative pressures indicate contraction due to attractive forces.</figcaption></figure>

MMT containing 0.67 e per unit cell. Ca²⁺ ions displayed lower diffusion compared to Na⁺ ions, a consistent trend attributed to the more stable coordination shell of Ca²⁺ as reported in prior research [30].

Interestingly, the diffusion pattern of Cu²⁺ ions is lower than that of Ca²⁺ ions in all directions and for all three hydration layers. However, when compared to Na⁺ ions, the diffusion behavior of Cu²⁺ ions is more similar to that of Ca²⁺ ions. This suggests that, despite Cu²⁺ ions having lower mobility than Ca²⁺ ions, they share more comparable mobility characteristics with Ca²⁺ ions than with Na⁺ ions in the context of different hydration layers. The higher diffusion rate of Na⁺ ions is due to their single positive charge, which results in weaker interactions with water molecules compared to the divalent Ca²⁺ and Cu²⁺ ions. Ca²⁺ and Cu²⁺ have similar diffusion behaviors because of their similar charges, but the slightly lower diffusion rate of Cu²⁺ can be attributed to its somewhat smaller hydrated ionic radius, which leads to stronger interactions with the surrounding water molecules, thus reducing its mobility.

## 4. Discussions and conclusions

In this work, we employed DFT calculations for Cu-MMT and extended the ClayFF force field to include Cu²⁺ interactions, allowing for

the accurate modeling of Cu-MMT systems. We examined the interaction energies between two MMT with three different cations including Na⁺, Ca²⁺, and Cu²⁺. Our calculations suggest that, in pure water solution, the swelling pressure induced by Cu-MMT falls between that of Ca-MMT and Na-MMT, a significant observation with implications for DGR applications. Our study suggests that exerted swelling pressure caused by Cu-MMT on the container is comparable to those of Ca- and Na-MMT. Also, Cu-MMT exhibits substantial swelling capabilities similar to Na- and Ca-MMT, effectively hindering bacterial activities, which is crucial for the protection of the container and therefore the containment of used nuclear fuel. Of note, in a DGR system, the actual conditions might involve the presence of saline solutions, which could influence the Cu ion charge state, the swelling pressure induced by Cu-MMT and its implications for DGR applications.

Furthermore, our MSD studies demonstrated that Na⁺ ions exhibit higher mobility than Ca²⁺ and Cu²⁺ ions, particularly in the XY plane. We observed that the diffusion rate increases with the hydration layer for all ions and in all directions. This observation indicates that additional water layers facilitate greater ion mobility, likely due to reduced interactions with the clay surface as hydration increases.

We found that the diffusion of Cu²⁺ ions was lower than that of Ca²⁺ ions in all directions and for all three hydration layers. This behavior is expected based on their different chemical properties. Our calculations
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Mean square displacement (MSD) calculations for interlayer Na⁺, Ca²⁺, and Cu²⁺ ions in one, two, and three hydrated layers. Columns represent the number of hydrated layers: first column (1 layer), second column (2 layers), and third column (3 layers). Rows display total MSD, as well as x, y, and z components of the MSD.</figcaption></figure>

<caption>**Table 3**</caption>

Comparison of self-diffusion coefficients (in $10^{-10}$ m²/s) for Cu²⁺, Ca²⁺, and Na⁺ at three different hydration levels in the interlayer regions.

<table><thead><tr><th>Layers of hydration</th><th>D<sub>Na+</sub></th><th>D<sub>Ca²+</sub></th><th>D<sub>Cu²+</sub></th></tr></thead><tbody><tr><td>1-layer hydrated</td><td></td><td></td><td></td></tr><tr><td>X-direction</td><td>0.99 ± 6.25 × 10<sup>-5</sup></td><td>0.09 ± 1.94 × 10<sup>-5</sup></td><td>1.61 × 10<sup>-2</sup> ± 1.62 × 10<sup>-6</sup></td></tr><tr><td>Y-direction</td><td>1.09 ± 6.57 × 10<sup>-5</sup></td><td>0.16 ± 3.39 × 10<sup>-5</sup></td><td>1.11 × 10<sup>-2</sup> ± 6.05 × 10<sup>-6</sup></td></tr><tr><td>Total</td><td>1.04 ± 3.57 × 10<sup>-5</sup></td><td>0.13 ± 2.19 × 10<sup>-5</sup></td><td>1.32 × 10<sup>-2</sup> ± 3.36 × 10<sup>-6</sup></td></tr><tr><td>2-layer hydrated</td><td></td><td></td><td></td></tr><tr><td>X-direction</td><td>4.98 ± 4.56 × 10<sup>-4</sup></td><td>1.98 ± 3.50 × 10<sup>-4</sup></td><td>0.53 ± 5.47 × 10<sup>-5</sup></td></tr><tr><td>Y-direction</td><td>3.90 ± 2.68 × 10<sup>-4</sup></td><td>1.38 ± 5.14 × 10<sup>-4</sup></td><td>0.42 ± 7.98 × 10<sup>-5</sup></td></tr><tr><td>Total</td><td>4.44 ± 3.2 × 10<sup>-4</sup></td><td>1.58 ± 1.2 × 10<sup>-4</sup></td><td>0.49 ± 6.93 × 10<sup>-5</sup></td></tr><tr><td>3-layer hydrated</td><td></td><td></td><td></td></tr><tr><td>X-direction</td><td>7.34 ± 3.34 × 10<sup>-4</sup></td><td>2.31 ± 4.24 × 10<sup>-4</sup></td><td>1.15 ± 2.33 × 10<sup>-4</sup></td></tr><tr><td>Y-direction</td><td>7.74 ± 5.64 × 10<sup>-4</sup></td><td>2.32 ± 3.96 × 10<sup>-4</sup></td><td>2.25 ± 4.15 × 10<sup>-4</sup></td></tr><tr><td>Total</td><td>7.53 ± 1.83 × 10<sup>-4</sup></td><td>2.29 ± 3.79 × 10<sup>-4</sup></td><td>1.72 ± 1.42 × 10<sup>-4</sup></td></tr></tbody></table>

showed that, in the event of copper corrosion from the container, these Cu²⁺ ions can be captured by the surrounding bentonite.

In this research, we considered Cu²⁺ as the charge state of copper ions in the clay systems. However, copper may exist in other charge states. This assumption is empowered by the results of energy profile analysis of more than 40 L-J potentials for Cu²⁺ along the perpendicular direction of MMT, which two of the best results are presented in Fig. 2. Despite testing a variety of potentials, we found that our L-J potential captures the well depth and convexity of the curvature around the equilibrium better than other tested L-J potentials, but it still does not perfectly match the DFT results. This discrepancy suggests that the charge of the copper ion may not be exactly 2⁺ in the presence of dry MMT. Furthermore, it is possible that the charge of the copper ion varies with its distance from MMT. In future, we will explore the charge

distribution of copper ion using Bader charge analysis in dry and wet clay systems to gain deeper insight into these phenomena.

Overall, our findings provide a better understanding of the behavior of interlayer cations in clay systems. The extension of the ClayFF force field to include Cu²⁺ interactions enhances its versatility and applicability in modeling a broader range of clay systems. In the future, we aim to use the PMFs in this study to develop coarse-grained mesoscales of hydrated Cu-MMT, Ca-MMT and Na-MMT and examine their mechanical properties. An interesting extension of this work would involve including a PMF calculation where the rotation of one platelet around the Z-axis is relaxed, allowing for the actual minimum PMF to be obtained. This could potentially alter the observed swelling pressures and provide a more accurate representation of the interaction characteristics.
CRediT authorship contribution statement

**Yalda Pedram:** Conceptualization, Data curation, Formal analysis, Investigation, Validation, Visualization, Writing – original draft. **Yaoting Zhang:** Conceptualization, Formal analysis, Methodology. **Scott Briggs:** Conceptualization, Writing – review & editing. **Chang Seok Kim:** Validation, Writing – review & editing. **Laurent Brochard:** Formal analysis, Methodology, Validation, Writing – review & editing. **Andrey G. Kalinichev:** Conceptualization, Methodology. **Laurent Karim Béland:** Conceptualization, Formal analysis, Funding acquisition, Project administration, Supervision, Validation, Writing – review & editing.

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

This work is financially supported by the Nuclear Waste Management Organization, Canada, the Natural Science and Engineering Research Council (NSERC), Canada, and Mitacs, Canada. We also thank the Digital Research Alliance of Canada and Centre for Advanced Computing at Queen's University for generous allocation of computer resources.

Data availability

Data will be made available on request.

References

[1] S. Guggenheim, J.M. Adams, D. Bain, F. Bergaya, M.F. Brigatti, V. Drits, M.L. Formoso, E. Galán, T. Kogure, H. Stanjek, Summary of recommendations of nomenclature committees relevant to clay mineralogy: report of the Association Internationale pour l'Etude des Argiles (AIPEA) Nomenclature Committee for 2006, Clays Clay Miner. 54 (2006) 761-772.

[2] H.H. Murray, Bentonite applications, Dev. Clay Sci. 2 (2006) 111-130.

[3] A.G. Clem, R.W. Doehler, Industrial applications of bentonite, Clays Clay Miner. 10 (1961) 272-283.

[4] P. Abootalebi, G. Siemens, Thermal properties of engineered barriers for a Canadian deep geological repository, Can. Geotech. J. 55 (6) (2018) 759-776.

[5] F. Uddin, Clays, nanoclays, and montmorillonite minerals, Metall. Mater. Trans. A 39 (12) (2008) 2804-2814.

[6] A.T. Bourgoyne, K.K. Millheim, M.E. Chenevert, F.S. Young, et al., Applied Drilling Engineering, vol. 2, Society of Petroleum Engineers Richardson, 1986.

[7] M.G. Temraz, I. Hassanien, Mineralogy and rheological properties of some Egyptian bentonite for drilling fluids, J. Nat. Gas Sci. Eng. 31 (2016) 791-799.

[8] R.O. Afolabi, O.D. Orodu, V.E. Efeobvokhan, Properties and application of Nigerian bentonite clay deposits for drilling mud formulation: recent advances and future prospects, Appl. Clay Sci. 143 (2017) 39-49.

[9] D.D. Carr, A.F. Alsobrook, G. Austin, Industrial Minerals and Rocks, Society for mining metallurgy and exploration, 1994.

[10] J.-H. Park, H.-J. Shin, M.H. Kim, J.-S. Kim, N. Kang, J.-Y. Lee, K.-T. Kim, J.I. Lee, D.-D. Kim, Application of montmorillonite in bentonite as a pharmaceutical excipient in drug delivery systems, J. Pharm. Investig. 46 (2016) 363-375.

[11] R.H.S. Robertson, Fullers earth: a history of calcium montmonillonite, 1986.

[12] E. Alvarez-Ayuso, A. García-Sánchez, Removal of heavy metals from waste waters by natural and na-exchanged bentonites, Clays Clay Miner. 51 (5) (2003) 475-480.

[13] X. Liu, C. Tournassat, S. Grangeon, A.G. Kalinichev, Y. Takahashi, M. Marques Fernandes, Molecular-level understanding of metal ion retention in clay-rich materials, Nat. Rev. Earth Environ. 3 (7) (2022) 461-476.

[14] P. Sellin, O.X. Leupin, The use of clay as an engineered barrier in radioactive-waste management-a review, Clays Clay Miner. 61 (6) (2013) 477-498.

[15] B. Grambow, Geological disposal of radioactive waste in clay, Elements 12 (4) (2016) 239-245.

[16] N. Giroud, Y. Tomonaga, P. Wersin, S. Briggs, F. King, T. Vogt, N. Diomidis, On the fate of oxygen in a spent fuel emplacement drift in Opalinus Clay, Appl. Geochem. 97 (2018) 270-278.

[17] SKB, Post-Closure Safety for the Final Repository for Spent Nuclear Fuel at Farsmark, Buffer, Backfill and Closure Process Report, PSAR Version, 2022.

[18] J. Aromaa, M. Kekkonen, M. Mousapour, A. Jokilaakso, M. Lundström, The oxidation of copper in air at temperatures up to 100 C, Corros. Mater. Degrad. 2 (4) (2021) 625-640.

[19] F. Uddin, et al., Montmorillonite: An Introduction to Properties and Utilization, vol. 817, IntechOpen London, UK, 2018.

[20] J.R. Booker, R. Brachman, R. Quigley, R.K. Rowe, Barrier Systems for Waste Disposal Facilities, Crc Press, 2004.

[21] I. Odom, Smectite clay minerals: properties and uses, Phil. Trans. R. Soc. A 311 (1517) (1984) 391-409.

[22] R. Anderson, I. Ratcliffe, H. Greenwell, P. Williams, S. Cliffe, P. Coveney, Clay swelling—a challenge in the oilfield, Earth-Sci. Rev. 98 (3-4) (2010) 201-216.

[23] R.T. Cygan, J.A. Greathouse, A.G. Kalinichev, Advances in clayf molecular simulation of layered and nanoporous materials and their aqueous interfaces, J. Phys. Chem. C 125 (32) (2021) 17573-17589.

[24] K. Norrish, Crystalline swelling of montmorillonite: manner of swelling of montmorillonite, Nature 173 (4397) (1954) 256-257.

[25] K. Norrish, The swelling of montmorillonite, Discuss. Faraday Soc. 18 (1954) 120-134.

[26] D. Fink, G. Thomas, X-ray studies of crystalline swelling in montmorillonites, Soil Sci. Am. J. 28 (6) (1964) 747-750.

[27] Y. Zhang, J. Oestreicher, W.J. Binns, S. Briggs, C.S. Kim, L.K. Béland, A coarse-grained interaction model for sodium dominant montmorillonite, Langmuir 38 (43) (2022) 13226-13237.

[28] D. Dixon, Review of the T-h-m-c properties of MX-80 bentonite, Nuclear Waste Management Organization, NWMO TR-2019-07, Toronto, 2019.

[29] I. Ravina, P.F. Low, Change of b-dimension with swelling of montmorillonite, Clays Clay Miner. 25 (1977) 201-204.

[30] L. Zhang, X. Lu, X. Liu, J. Zhou, H. Zhou, Hydration and mobility of interlayer ions of (Na X, Ca Y)-montmorillonite: a molecular dynamics study, J. Phys. Chem. C 118 (51) (2014) 29811-29821.

[31] W.R. Foster, J. Savins, J. Waite, Lattice expansion and rheological behavior relationships in water-montmorillonite systems, Clays Clay Miner. 3 (1954) 296-316.

[32] A. Meleshyn, C. Bunnenberg, The gap between crystalline and osmotic swelling of Na-montmorillonite: A Monte Carlo study, J. Chem. Phys. 122 (3) (2005).

[33] B. Derjaguin, L. Landau, Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes, Prog. Surf. Sci. 43 (1-4) (1993) 30-59.

[34] R.M. Pashley, J.N. Israelachvili, Molecular layering of water in thin films between mica surfaces and its relation to hydration forces, J. Colloid Interface Sci. 101 (2) (1984) 511-523.

[35] J.N. Israelachvili, R.M. Pashley, Molecular layering of water at surfaces and origin of repulsive hydration forces, Nature 306 (5940) (1983) 249-250.

[36] R. Jellander, S. Marčelja, J. Quirk, Attractive double-layer interactions between calcium clay particles, J. Colloid Interface Sci. 126 (1) (1988) 194-211.

[37] J. Du, A. Zhou, X. Lin, Y. Bu, J. Kodikara, Revealing expansion mechanism of cement-stabilized expansive soil with different interlayer cations through molecular dynamics simulations, J. Phys. Chem. C 124 (27) (2020) 14672-14684.

[38] R.T. Cygan, J.-J. Liang, A.G. 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.

[39] T.A. Ho, L.J. Criscenti, J.A. Greathouse, Revealing transition states during the hydration of clay minerals, J. Phys. Chem. Lett. 10 (13) (2019) 3704-3709.

[40] T. Honorio, L. Brochard, M. Vandamme, Hydration phase diagram of clay particles from molecular simulations, Langmuir 33 (44) (2017) 12766-12776.

[41] L. Sun, J.T. Hirvi, T. Schatz, S. Kasa, T.A. Pakkanen, Estimation of montmorillonite swelling pressure: A molecular dynamics approach, J. Phys. Chem. C 119 (34) (2015) 19863-19868.

[42] B. Akinwunmi, J.T. Hirvi, S. Kasa, T.A. Pakkanen, Swelling pressure of Na- and Ca-montmorillonites in saline environments: A molecular dynamics study, Chem. Phys. 528 (2020) 110511.

[43] B. Akinwunmi, Molecular Dynamics Study of the Bentonite Barrier: Effects of Salinity and Temperature on the Swelling Pressure (Ph.D. thesis), Itä-Suomen yliopisto, 2020.

[44] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, et al., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter. 21 (39) (2009) 395502.

[45] P.E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50 (24) (1994) 17953.

[46] A. Dal Corso, Pseudopotentials periodic table: From h to pu, Comput. Mater. Sci. 95 (2014) 337-350.

[47] We used the pseudopotentials Mg.pbe-spnl-kjpaw-psl.1.0.0UPF, Al.pbe-nl-kjpaw-psl.1.0.0UPF, O.pbe-n-kjpaw-psl.1.0.0UPF, H.pbe-kjpaw-psl.1.0.0UPF, Si.pbe-nl-kjpaw-psl.1.0.0UPF, Cu.pbe-dn-kjpaw-psl.1.0.0UPF, and Ca.pbe-spn-kjpaw-psl.1.0.0UPF from (http://www.quantum-espresso.org).
[48] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (18) (1996) 3865.

[49] 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 (15) (2010).

[50] M. Holmboe, Atom: A MATLAB package for manipulation of molecular systems, Clays Clay Miner. 67 (5) (2019) 419–426.

[51] H.J. Monkhorst, J.D. Pack, Special points for brillouin-zone integrations, Phys. Rev. B 13 (12) (1976) 5188.

[52] R. Fletcher, Practical Methods of Optimization, John Wiley & Sons, 2000.

[53] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017.

[54] H.J. Berendsen, J.P. Postma, W.F. van Gunsteren, J. Hermans, Interaction models for water in relation to protein hydration, in: Intermolecular Forces: Proceedings of the Fourteenth Jerusalem Symposium on Quantum Chemistry and Biochemistry, Jerusalem, Israel, April 13–16, 1981, Springer, 1981, pp. 331–342.

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

[56] A. Dufresne, J. Arayro, T. Zhou, K. Ioannidou, F.-J. Ulm, R. Pellenq, L.K. Béland, Atomistic and mesoscale simulation of sodium and potassium adsorption in cement paste, J. Chem. Phys. 149 (7) (2018).

[57] J.A. Greathouse, J.S. Durkin, J.P. Larentzos, R.T. Cygan, Implementation of a morse potential to model hydroxyl behavior in phyllosilicates, J. Chem. Phys. 130 (13) (2009).

[58] M. Povreau, J.A. Greathouse, R.T. Cygan, A.G. 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.

[59] L.N. Lammers, I.C. Bourg, M. Okumura, K. Kolluri, G. Sposito, M. Machida, Molecular dynamics simulations of cesium adsorption on illite nanoparticles, J. Colloid Interface Sci. 490 (2017) 608–620.

[60] M. Povreau, J.A. Greathouse, R.T. Cygan, A.G. 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.

[61] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1) (1995) 1–19.

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

[63] M.R. Hestenes, E. Stiefel, et al., Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 49 (6) (1952) 409–436.

[64] S. Majidi, H. Erfan-Niya, J. Azamat, E.R. Cruz-Chú, J.H. Walther, The performance of a C2N membrane for heavy metal ions removal from water under external electric field, Sep. Purif. Technol. 289 (2022) 120770.

[65] V. Moses, Ö.T. Bishop, K.A. Lobb, The evaluation and validation of copper (II) force field parameters of the Auxiliary Activity family 9 enzymes, Chem. Phys. Lett. 678 (2017) 91–97.

[66] K. Anitha, S. Namsani, J.K. Singh, Removal of heavy metal ions using a functionalized single-walled carbon nanotube: a molecular dynamics study, J. Phys. Chem. A 119 (30) (2015) 8349–8358.

[67] P. Li, K.M. Merz Jr., Taking into account the ion-induced dipole interaction in the nonbonded model of ions, J. Chem. Theory Comput. 10 (1) (2014) 289–297.

[68] P. Guan, D. Mckenzie, B. Pailthorpe, MD simulations of ag film growth using the lennard-jones potential, J. Phys.: Condens. Matter. 8 (45) (1996) 8753.

[69] K. Louw, B. Bradshaw-Hajek, J. Hill, Interaction of ferric ions with europium metal organic framework and application to mineral processing sensing, Phil. Mag. 102 (13) (2022) 1231–1246.

[70] T. Akiner, H. Ertürk, K. Atalik, Prediction of thermal conductivity and shear viscosity of water-cu nanofluids using equilibrium molecular dynamics, in: ASME International Mechanical Engineering Congress and Exposition, vol. 56369, American Society of Mechanical Engineers, 2013, p. V08CT09A012.

[71] H. Khanmohammadi, B. Bayati, J. Rahbar-Shahrouzi, A.-A. Babaluo, A. Ghorbani, Molecular simulation of the ion exchange behavior of Cu2+, Cd2+ and Pb2+ ions on different zeolites exchanged with sodium, J. Environ. Chem. Eng. 7 (3) (2019) 103040.

[72] J. Torras, C. Aleman, Determination of new Cu+, Cu2+, and Zn2+ lennard-jones ion parameters in acetonitrile, J. Phys. Chem. B 117 (36) (2013) 10513–10522.

[73] P. Li, B.P. Roberts, D.K. Chakravorty, K.M. Merz Jr., Rational design of particle mesh Ewald compatible lennard-jones parameters for 2 metal cations in explicit solvent, J. Chem. Theory Comput. 9 (6) (2013) 2733–2748.

[74] Q. Liao, S.C.L. Kamerlin, B. Strodel, Development and application of a non-bonded Cu2+ model that includes the Jahn–Teller effect, J. Phys. Chem. Lett. 6 (13) (2015) 2657–2662.

[75] E. Duffour, P. Malfreyt, MD simulations of the collision between a copper ion and a polyethylene surface: an application to the plasma–insulating material interaction, Polymer 45 (13) (2004) 4565–4575.

[76] C. Zhang, L. Liu, H. Jia, Molecular dynamics study of the exchange processes of heavy metals into montmorillonite: Characterization of hydrated edge surfaces and dynamic exchange mechanism, Appl. Geochim. 150 (2023) 105587.

[77] A.W. Knight, P. Ilani-Kashkouli, J.A. Harvey, J.A. Greathouse, T.A. Ho, N. Kabengi, A.G. Ilgen, Interfacial reactions of Cu (II) adsorption and hydrolysis driven by nano-scale confinement, Environ. Sci. Nano 7 (1) (2020) 68–80.

[78] S. Izvekov, G.A. Voth, Car–Parrinello molecular dynamics simulation of liquid water: New results, J. Chem. Phys. 116 (23) (2002) 10372–10376.

[79] D. Trzesniak, A.-P.E. Kunz, W.F. van Gunsteren, A comparison of methods to compute the potential of mean force, ChemPhysChem 8 (1) (2007) 162–169.

[80] J. Kästner, Umbrella sampling, Wiley Interdiscip. Rev. Comput. Mol. Sci. 1 (6) (2011) 932–942.

[81] M. Sprik, G. Ciccotti, Free energy from constrained molecular dynamics, J. Chem. Phys. 109 (18) (1998) 7737–7744.

[82] V.K. Voora, W. Al-Saidi, K.D. Jordan, Density functional theory study of pyrophyllite and M-montmorillonites (M=Li, Na, K, Mg, and Ca): Role of dispersion interactions, J. Phys. Chem. A 115 (34) (2011) 9695–9703.

[83] G.W. Marini, K.R. Liedl, B.M. Rode, Investigation of Cu2+ hydration and the Jahn–Teller effect in solution by QM/MM Monte Carlo simulations, J. Phys. Chem. A 103 (51) (1999) 11387–11393.

[84] D.T. Bowron, M. Amboage, R. Boada, A. Freeman, S. Hayama, S. Díaz-Moreno, The hydration structure of Cu 2+: more tetrahedral than octahedral? RSC Adv. 3 (39) (2013) 17803–17812.

[85] C.F. Schwenk, B.M. Rode, Extended ab initio quantum mechanical/molecular mechanical molecular dynamics simulations of hydrated Cu 2+, J. Chem. Phys. 119 (18) (2003) 9523–9531.

[86] A.C. Van Duin, V.S. Bryantsev, M.S. Diallo, W.A. Goddard, O. Rahaman, D.J. Doren, D. Raymand, K. Hermansson, Development and validation of a ReaFF reactive force field for Cu cation/water interactions and copper metal/metal oxide/metal hydroxide condensed phases, J. Phys. Chem. A 114 (35) (2010) 9507–9514.

[87] H. Ohtaki, T. Radnai, Structure and dynamics of hydrated ions, Chem. Rev. 93 (3) (1993) 1157–1204.

[88] S. Amira, D. Spångberg, K. Hermansson, Distorted five-fold coordination of Cu 2+(aq) from a Car–Parrinello molecular dynamics simulation, Phys. Chem. Chem. Phys. 7 (15) (2005) 2874–2880.

[89] M. Magini, Coordination of copper (II). Evidence of the Jahn–Teller effect in aqueous perchlorate solutions, Inorg. Chem. 21 (4) (1982) 1535–1538.

[90] A. Musinu, G. Paschina, G. Piccaluga, M. Magini, Coordination of copper (II) in aqueous copper sulfate solution, Inorg. Chem. 22 (8) (1983) 1184–1187.

[91] A. Pasquarello, I. Petri, P.S. Salmon, O. Parisel, R. Car, É. Tóth, D.H. Powell, H.E. Fischer, L. Helm, A.E. Merbach, First solvation shell of the Cu (II) aqua ion: evidence for fivefold coordination, Science 291 (5505) (2001) 856–859.

[92] B. Beagley, A. Eriksson, J. Lindgren, I. Persson, L. Pettersson, M. Sandstrom, U. Wahlgren, E. White, A computational and experimental study on the Jahn–Teller effect in the hydrated copper (II) ion. Comparisons with hydrated nickel (II) ions in aqueous solution and solid Tutton’s salts, J. Phys.: Condens. Matter. 1 (13) (1989) 2395.

[93] G. Clark, R. Grim, W. Bradley, A study of the behavior of montmorillonite upon wetting, Zeitschrift Kristallographie–Cryst. Mater. 97 (1–6) (1937) 216–222.

[94] M. Antognozzi, A. Humphris, M. Miles, Observation of molecular layering in a confined water film and study of the layers viscoelastic properties, Appl. Phys. Lett. 78 (3) (2001) 300–302.

[95] L. Brochard, T. Honório, M. Vandamme, M. Bornert, M. Peigney, Nanoscale origin of the thermo-mechanical behavior of clays, Acta Geotech. 12 (2017) 1261–1279.

[96] L. Brochard, Swelling of montmorillonite from molecular simulations: Hydration diagram and confined water properties, J. Phys. Chem. C 125 (28) (2021) 15527–15543.

[97] B. Akinwunmi, L. Sun, J.T. Hirvi, S. Kasa, T.A. Pakkanen, Influence of temperature on the swelling pressure of bentonite clay, Chem. Phys. 516 (2019) 177–181.

[98] S.S. Agus, T. Schanz, A method for predicting swelling pressure of compacted bentonites, Acta Geotech. 3 (2) (2008) 125–137.

[99] S. Kumpulainen, L. Kiviranta, Mineralogical, Chemical and Physical Study of Potential Buffer and Backfill Materials from ABM. Test Package 1, Tech. Rep., Posiva Oy, 2011.

[100] O. Karnland, S. Olsson, U. Nilsson, Mineralogy and Sealing Properties of Various Bentonites and Smectite-Rich Clay Materials, Tech. Rep., Swedish Nuclear Fuel and Waste Management Co., 2006.

[101] F. Bucher, U. Spiegel, Quelldruck von hochverdichteten Bentoniten, Nagra, 1984.

[102] P. Sitz, Materialuntersuchungen für Mehrkomponentensysteme auf Ton/Bentonit Basis für Dichtung und Lastabtrag, mit hohem Ruckhaltvermögen, für den langzeitsicheren Verschlußvon UTD und Lastabtrag, 1997.

[103] M. Holmboe, I.C. Bourg, Molecular dynamics simulations of water and sodium diffusion in smectite interlayer nanopores as a function of pore size and temperature, J. Phys. Chem. C 118 (2) (2014) 1001–1013.

[104] J.A. Greathouse, R.T. Cygan, J.T. Friedrich, G.R. Jerauld, Molecular dynamics simulation of diffusion and electrical conductivity in montmorillonite interlayers, J. Phys. Chem. C 120 (3) (2016) 1640–1649.