# Atrazine desorption mechanism from an hydrated calcium montmorillonite-A DFT molecular dynamics study

Quentin Desdion, Fabienne Bessac, Sophie Hoyau

► To cite this version:

Quentin Desdion, Fabienne Bessac, Sophie Hoyau. Atrazine desorption mechanism from an hydrated calcium montmorillonite-A DFT molecular dynamics study. International Journal of Molecular Sciences, 2024, 25 (3), pp.1604. 10.3390/ijms25031604. hal-04502975

HAL Id: hal-04502975

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

Submitted on 2 Jun 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

# Atrazine Desorption Mechanism from an Hydrated Calcium Montmorillonite—A DFT Molecular Dynamics Study

Quentin Desdion <fnref n="1" />,<fnref n="2" />, Fabienne Bessac <fnref n="1" />,<fnref n="2" />,* and Sophie Hoyau <fnref n="1" />,†

**Abstract:** Atrazine is one of the most widely used herbicide molecules in the triazine family. Despite its interdiction in the European Union in 2004, atrazine and its main degradation products remain among the most frequently found molecules in freshwater reservoirs in many European Union countries. Our study aims in obtaining insight into the desorption process of atrazine from the main soil absorbent material: clay. Constrained Molecular Dynamics simulations within the Density Functional Theory framework allow us to obtain a free energy desorption profile of atrazine from a Ca<fnref n="2" />⁺-montmorillonite surface. The results are interpreted in terms of atrazine inclination to the clay surface and moreover, in terms of hydration states of the cations present in the clay interlayer as well as the hydration state of the atrazine. The desorption mechanism is driven by atrazine alkyl groups and their sizes because of dispersion stabilizing effects. The highest barrier corresponds to the loss of the isopropyl interaction with the surface.

**Keywords:** atrazine; clay; desorption; hydration; DFT; molecular dynamics

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

check for
updates

**Citation:** Desdion, Q.; Bessac, F.; Hoyau, S. Atrazine Desorption Mechanism from an Hydrated Calcium Montmorillonite—A DFT Molecular Dynamics Study. *Int. J. Mol. Sci.* **2024**, *25*, 1604. https://doi.org/10.3390/ijms25031604

Academic Editor: Andrzej Grzybowski

Received: 21 December 2023

Revised: 18 January 2024

Accepted: 22 January 2024

Published: 27 January 2024

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

**Copyright:** © 2024 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

The sorption of pesticides and, more broadly, of organic contaminants, is an important factor to study in soil and water contamination processes, as well as in the development of methods for removing pollutant molecules from these environmental compartments. A wide range of studies concerns the investigation of adsorption and desorption of pesticides from soils or absorbent materials such as biochar, clays or activated carbon [1–3]. Among the many pesticide molecules, atrazine (see Figure 1) is the active ingredient in herbicides.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Structural formula of the atrazine herbicide, C<sub>8</sub>H<sub>14</sub>ClN<sub>5</sub>, 2-isopropylamino-4-ethylamino-6-chloro-1,3,5-triazine, and numbering of the atoms. Ce is the label for the ethyl carbon atom linked to the amino group and Ci is the label for the isopropyl carbon atom linked to the amino group.</figcaption></figure>

It is one of the most widely used herbicide molecules in the triazine family. Atrazine is a selective herbicide that controls the development of broadleaf weeds and grasses in major

<footnote n="1">¹ Laboratoire de Chimie et Physique Quantiques, Université de Toulouse, Université Paul Sabatier, Toulouse III, CNRS (UMR 5626), 118 Route de Narbonne, F-31062 Toulouse, France; qdesdion@irsamc.ups-tlse.fr (Q.D.); sophie.hoyau@irsamc.ups-tlse.fr (S.H.)</footnote>

<footnote n="2">² Ecole d'Ingénieurs de Purpan, Université de Toulouse, Toulouse INP, 75 Voie du TOEC, BP 57611, Cedex 03, F-31076 Toulouse, France</footnote>

<footnote>* Correspondence: fabienne.bessac@irsamc.ups-tlse.fr; Tel.: +33-762-062-786</footnote>
crops such as maize, soybean, citrus fruits and sugarcane. Atrazine and its main degradation products are of concern due to their potential endocrine and carcinogenic activity [4]. The introduction of atrazine dates back to the late 1950s. Since its introduction, atrazine has heavily been used. Thus, atrazine and its metabolites are widely present in the environment, contaminating ground and surface water, entering the atmosphere and suspected of persisting for decades in soils. In Europe, the use of atrazine was curtailed in the early 1990s due to excessively high concentrations found in drinking water [5]. In 2004, the European Union finally banned its use. However, atrazine remains one of the most frequently found molecules in freshwater reservoirs in many European Union countries, including France and Germany. Clays in soils play an important role in immobilizing many organic contaminants, such as atrazine and its metabolites [6]. Thus, atrazine has extensively been studied using experimental and theoretical resources and techniques. For instance, G. Abate et al. conducted adsorption and desorption experiments. They described the modification of the clay minerals, vermiculite and montmorillonite, by intercalating Fe(III) polymers with the aim of removing atrazine and its metabolites: deethylatrazine, deisopropylatrazine, and hydroxyatrazine from aqueous solution. The results of the desorption studies in water allow for the depiction of the following affinity order for K⁺ homoionic montmorillonite clay: hydroxyatrazine > deisopropylatrazine ≈ atrazine ≈ deethylatrazine [7]. In 2011, N. D. Jablonowski et al. wrote a short commentary article to alert about the use of atrazine regarding its persistence and toxicity and those of its metabolites [4]. A wide variety of theoretical and experimental studies have been carried out on atrazine, some of which combine theory and experience. We will comment here on some of the most recently published. Y. Han et al.'s computer simulation results aim for serving as a guideline for experimental design in a study to detect the trace level of atrazine residue in the environment and crops [8]. In another study similar to that of Y. Cheng et al., DFT calculations provided an interesting tool to reinforce the experimental results. Here, electrostatic potential showed evidence that the good adsorption performance of nitrogen (N)-doped cellulose biochar for atrazine mainly depends on chemisorption, and π – π electron donor-acceptor interaction [2]. In a work carried out by the V. N. Freire group, the binding of atrazine to human serum albumin was characterized. Describing those interactions is clinically relevant for future pharmacological and toxicological studies. The binding was investigated by using fluorescence spectroscopy and performing simulations using molecular docking, classical molecular dynamics and quantum biochemistry based on density functional theory (DFT) [9]. Similarly, Y. Vieira et al. explored adsorption mechanisms of two pesticides, atrazine and diuron, onto a model of activated carbon and used conceptual DFT descriptors to analyze, confirm or complete their experimental results. The atrazine molecule changes its interaction orientation from inclined to parallel to the activated carbon surface, resulting in higher electrostatic repulsion, reflecting on its low increase in adsorption capacity as the temperature increases [3]. In 2021, D. Tunega and coworkers worked on the stability of atrazine–smectite intercalates using DFT simulations and implementing experimental data. Smectites of this study were montmorillonite and beidellite. The calculations revealed an arrangement of the atrazine in the interlayer space parallel to the surface of both smectites. The stability of beidellite structures was higher. Moreover, the authors showed that organically modified smectite increased the fixation of atrazine [10].

Our study is focused on the desorption of atrazine from a Ca²⁺-montmorillonite surface (Mont) in the presence of water. Previously, we published the results of simulations on the desorption of a fungicide, the fenhexamid, from the same hydrated surface [11]. In this work, we first present the results: we will describe the free energy desorption profile before we discuss the free energy changes in terms of atrazine inclination to the Mont surface, dispersion effects and hydration state of both the cations and of the atrazine molecule. In the Materials and Methods part, we will present the model and the methods we used to investigate the desorption of the atrazine from Mont.
**2. Results and Discussion**

**2.1. Description of the Desorption Profile**

In Figure 2, the free energy desorption profile of Atra from surf is depicted. Five singularities are noticeable and labeled P1, M1, P2, M2 and P3. The zones delimited by the red horizontal lines enclose the structures considered as belonging to the minima and to the maxima. P1, P2 and P3 are free energy minima separated by two barriers evaluated from the smoothed curve at 1.5 and 11.5 kcal·mol⁻¹, a low and a high barrier, respectively. Taking the minimum free energy point for a well and the maximum free energy point for a maximum, these barriers are worth 2 and 13 kcal·mol⁻¹, respectively. M1 and M2 are maxima along the reaction coordinate ζ (ζ<sub>M1</sub> = 3.88 Å and ζ<sub>M2</sub> = 5.35 Å). P1 and P2 are quasi isoenergetic: P1 and P2 are centered on ζ values of 3.41 Å and 4.24 Å, respectively. In addition, P3 centered on ζ = 5.56 Å is close to M2 in free energy (≈0.5 kcal·mol⁻¹). After P3, the free energy increases with ζ without anymore clear singularity. To analyze this part of the profile, we define three ζ ranges 0.4 Å wide centered on 6.0, 7.0 and 8.0 Å.

In the Supporting Information (SI), information is detailed confirming the liquid phase of water all along the desorption area (see Figure S2a). Moreover, water molecules are not directly in contact with the Mont surface and the hydration of the surface stays similar all along the desorption process (see Figure S2b).

The compensating cations, in our case Ca²⁺, are initially located close to surf (<0.5 Å) and in the center of rings formed by six Ob atoms, thus they lie above an O sheet isomorphic substitution. Ca1 is defined as the cation interacting with Atra. The starting interaction mode is bidentate and involves Cl, the chlorine atom linked to C6 carbon atom and N1, the nitrogen atom of the triazine cycle bonded to C6 and the closest to the isopropyl group (see Figure 1). This complexation site corresponds to the lowest-energy isomer in the gas phase for the Atra-Ca²⁺ complex at the PBE level (unpublished results).

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

**Figure 2.** Free energy profile of atrazine desorption from the montmorillonite surface as a function of ζ, the distance from the center of mass of atrazine to the surface containing the oxygen atoms. Ob. Every black dot represents a structure visited during the 19 MD trajectories. The continuous purple line gives the free energy profile of the desorption after smoothing. P1, P2 and P3 are the minima or energetic wells of this profile. M1 and M2 are the maxima. The zones delimited by the red horizontal lines enclose the structures considered as belonging to the minima and to the maxima. Numbers on the top of the graph are the ζ values for each singularity.
## 2.2. Interaction of Atrazine with the Ca1 Cation

To analyze the complexation site of Ca1 on the atrazine molecule, we will refer to Figure 3a,b. In Figure 3a, the normalized Ca1-N1 distance distribution shows that for P1, M1 and P2, the interaction involves N1, a nitrogen atom of the triazine cycle. From M2 and larger $\zeta$ values, the Ca1-N1 interaction is lost. Furthermore, in Figure 3b, the normalized Ca1-Cl distance distribution is presented. For P2, this figure evidence the interaction of Ca1 with Cl, the chlorine atom of Atra (Ca1-Cl $\approx$ 3.2 Å), implying a bidentate complexation site between Ca1 and Atra. For P1 and M1, the Ca1-Cl distance is longer ($\approx$3.7–3.8 Å) leading to a monodate interaction site. From M2 and larger $\zeta$ values, all Ca1-Atra interactions are lost: for Ca1-N1 as well as the Ca1-Cl. The difference between P1 and P2, which are energetically quasi degenerated, concerns the interaction site; monodentate for P1 and bidentate for P2.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** (a) Ca1-N1 and (b) Ca1-Cl normalized distance distributions for each of the singularities P1, M1, P2, M2 and P3 and also for the $\zeta$ values 6.0, 7.0 and 8.0 Å.</figcaption></figure>

## 2.3. Atrazine Inclination with Respect to the Surface

We define the angle $\tilde{XClN3} = \theta$ to follow the evolution of the atrazine triazine cycle inclination with respect to the surface (see N3 in Figure 1). The coordinates of the ghost atom X are $(x_{N3}; y_{N3}; z_{Cl})$. For $\theta$ values close to zero, the triazine cycle is almost parallel to the Mont surface. When $\theta$ is positive, $z_{N3} > z_{Cl}$, whereas when $\theta$ is negative, $z_{N3} < z_{Cl}$ (see Figure 4a). In Figure 4b, $\theta$ normalized distribution is given for each singularity of the profile and also for the $\zeta$ values 6.0, 7.0 and 8.0 Å. From P1 to M2, $\theta$ increases.

As we can see in Table 1, the mean value of $\theta$ changes from 7.1° in P1 to 48.7° in M2. However, in M2, the distribution presents two peaks: one around 30° and the other one around 60°. P3, the next singularity, shows a single peak which mean value (56.7°) is also close to 60° as in M2. For the $\zeta$ values around 6.0, 7.0 and 8.0 Å, the mean $\theta$ value decreases from 40.3° to 27.1°. Around 6.0 and 7.0 Å, the angle distributions present two peaks centered on 0°, 50°, and on 20°, 40°, respectively. Around 8.0 Å, the distribution is wide from 0 to 50°, with a mean value of 27.1°. In the beginning of the desorption process, the Atra triazine cycle is almost parallel to the surface. $\theta$ oscillates around 0°. When it is negative, the N3 nitrogen atom is closer to the surface than Cl. Early in the process, negative angles can be found only in the P1 and M1 singularities. The increase in the $\theta$ angle from P1 to P3 moves together with a loss of dispersion effects due to the detachment of the functionalized triazine cycle from the surface but a conservation of the Cl-surface dispersion interaction. Thus, the largest variation in $\theta$ is observed between P2 ($\bar{\theta} = 20.8°$) and M2 ($\bar{\theta} = 48.7°$), which corresponds to the highest free energy barrier in the profile, i.e., 11.5 kcal·mol⁻¹. Beyond P3, those dispersion effects are also lost and $\theta$ angle varies more. The evolution of the atrazine inclination with respect to the surface, described by $\theta$, is correlated to the variation of atrazine and Ca1 hydration states.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** (a) $\theta$ angle definition: top view, $\theta > 0$; bottom view, $\theta < 0$. (b) $\theta$ normalized distributions for each of the singularities P1, M1, P2, M2 and P3 and also for the $\zeta$ values 6.0, 7.0 and 8.0 Å. $\theta$ is given in degrees.</figcaption></figure>

<table><caption>**Table 1.** The mean number of water molecules involved in the first hydration spheres of Ca1 ($HS_{Ca1}$), Ca2 ($HS_{Ca2}$), Atrazine ($HS_{Atra}$) and Atrazine–Ca1 entity ($HS_{Atra-Ca1}$) for each of the singularities P1, M1, P2, M2 and P3 and also for around the $\zeta$ values 6.0, 7.0 and 8.0 Å. $HS_{tot}$ is the mean number of water molecules around Atrazine, Ca1 and Ca2. Mean Ca1–Ca2 distances are also reported in Å, as well as the mean $\theta$ angle, noted $\bar{\theta}$ and given in degrees. When $HS_{Ca1}$ is included into $HS_{Atra}$, it is written in bold.</caption><thead><tr><th></th><th>HS<sub>tot</sub></th><th>HS<sub>Atra-Ca1</sub></th><th>Ca1 &sub; Atra<sup>a</sup></th><th>HS<sub>Atra</sub></th><th>HS<sub>Ca1</sub></th><th>HS<sub>Ca2</sub></th><th>Ca1–Ca2 (Å)</th><th>&theta;&#773; (°)</th></tr></thead><tbody><tr><td>P1</td><td>26.4</td><td><strong>21.9</strong></td><td><strong>1.2</strong></td><td><strong>21.9</strong></td><td><strong>1.2</strong></td><td>4.5</td><td>9.10</td><td>7.1</td></tr><tr><td>M1</td><td>28.6</td><td><strong>24.1</strong></td><td><strong>0.8</strong></td><td><strong>24.1</strong></td><td><strong>0.8</strong></td><td>4.5</td><td>8.57</td><td>13.1</td></tr><tr><td>P2</td><td>29.2</td><td><strong>23.9</strong></td><td><strong>1.4</strong></td><td><strong>23.9</strong></td><td><strong>1.4</strong></td><td>5.2</td><td>9.74</td><td>20.8</td></tr><tr><td>M2</td><td>33.1</td><td>30.4</td><td>2.1</td><td>30.3</td><td>2.1</td><td>2.7</td><td>10.34</td><td>48.7</td></tr><tr><td>P3</td><td>34.9</td><td>31.9</td><td>2.8</td><td>31.7</td><td>3.0</td><td>3.0</td><td>10.29</td><td>56.7</td></tr><tr><td>6.0 Å</td><td>37.8</td><td>33.6</td><td>2.7</td><td>33.3</td><td>3.0</td><td>4.2</td><td>10.33</td><td>40.3</td></tr><tr><td>7.0 Å</td><td>41.6</td><td>37.7</td><td>1.9</td><td>36.6</td><td>3.0</td><td>3.9</td><td>9.70</td><td>36.5</td></tr><tr><td>8.0 Å</td><td>41.6</td><td>38.1</td><td>1.1</td><td>35.8</td><td>3.3</td><td>3.5</td><td>10.00</td><td>27.1</td></tr></tbody></table>

## 2.4. Dispersion Effects and Atrazine and Cation Hydration States

In this part, we will discuss the hydration states of the atrazine molecule and of the $Ca^{2+}$ cations. In Table 2, the minimum, maximum and mean numbers of water molecules involved in the first Hydration Sphere (HS) of Atra are given.

<table><caption>**Table 2.** The minimum (MIN), maximum (MAX) and mean (MEAN) numbers of water molecules involved in the first HS of Atrazine for each of the singularities P1, M1, P2, M2 and P3 and also around the $\zeta$ values 6.0, 7.0 and 8.0 Å.</caption><thead><tr><th></th><th>MIN</th><th>MAX</th><th>MEAN</th></tr></thead><tbody><tr><td>P1</td><td>14</td><td>29</td><td>21.9</td></tr><tr><td>M1</td><td>15</td><td>32</td><td>24.1</td></tr><tr><td>P2</td><td>16</td><td>30</td><td>23.9</td></tr><tr><td>M2</td><td>24</td><td>35</td><td>30.3</td></tr><tr><td>P3</td><td>24</td><td>38</td><td>31.7</td></tr><tr><td>6.0 Å</td><td>25</td><td>41</td><td>33.3</td></tr><tr><td>7.0 Å</td><td>30</td><td>44</td><td>36.6</td></tr><tr><td>8.0 Å</td><td>26</td><td>43</td><td>35.8</td></tr></tbody></table>

<footnote><sup>a</sup> Number of water molecules in common between $HS_{Atra}$ and $HS_{Ca1}$.</footnote>
From P1 to P2, the number of water molecules in HS<sub>Atra</sub> varies from 14 to 32, with a mean value of [21.9; 24.1]. Atra hydration state for these singularities is linked with small θ values (see θ̄ in Table 1). Moreover, the water molecule(s) that hydrate(s) Ca1 are included into HS<sub>Atra</sub>, thus the complex Atra-Ca1 has the same hydration state as Atra. In Figure 5, snapshots representative of each singularity are depicted. They were chosen to reflect the mean θ angle and mean Atra and cation hydration states as much as possible.

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

**Figure 5.** Snapshots taken in the five singularities P1, M1, P2, M2 and P3 of the free energy profile and their ζ values. For a sake of clarity, only the Ob oxygen atoms are represented for the Mont surface. For the same reason, only the water molecules involved in the first HS of both atrazine and Ca1 are represented. HS<sub>Ca1</sub>, HS<sub>Ca2</sub> and HS<sub>Atra</sub> the number of water molecules involved in the first HS of each entity are noticed. In red, the θ angle is given for the snapshot. In parentheses, the mean values of HS<sub>Atra</sub> and θ are also given. Atom color codes: H, orangy beige; C, gray; N, blue; O, red; Cl, green; Ca, magenta.

For P1, θ̄ is low (7.1°), thus, Atra is almost parallel to the surface. The chosen P1 snapshot shows an angle of 6.1°, Atra parallel to the surface bonds to Ca1 within a monodentate manner via N1 (Ca1 – N1 = 2.411 Å). The water molecules of the Atra hydration sphere, HS<sub>Atra</sub> (=20 in P1 snapshot), are distributed around the pesticide but not between the surface and Atra, conserving the Atra-surf dispersion effects. Both alkyl groups, ethyl and isopropyl, lie at average distances to the surface similar to the reaction coordinate ζ (see Figure 6). The Ca1 hydration sphere in P1 includes one to two molecules; one being predominant (see Figure 7). At the same time, HS<sub>Ca2</sub> varies from four to six water molecules; four being preferred (see Figure 8).

For M1, the situation is quite similar to P1 with θ̄ slightly larger (13.1°) and thus, with HS<sub>Atra</sub> containing on average two more water molecules. In contrast to P1, Ce-surf and Ci-surf do not follow the same tendency: Ci-surf remains close to ζ, whereas Ce-surf becomes longer reducing the dispersion effects with the surface but allowing the entrance of water molecules in the space thus created (see Figure 6). In the M1 snapshot, we can observe that the monodentate interaction site between Atra and Ca1 is maintained (Ca1 – N1 = 2.683 Å). Concerning HS<sub>Ca1</sub> and HS<sub>Ca2</sub>, more hydration states are accessible than in P1: zero to two, and two to six, respectively. The most frequent combination is HS<sub>Ca1</sub> = 1 and HS<sub>Ca2</sub> = 5 (see Figures 7 and 8).

In P2, the θ angle keeps on increasing (θ̄ = 20.8°). Ci-surf distance distribution stays similar to M1. The dispersion interaction between the largest alkyl group, i.e., isopropyl, and the surface stands still. On the other hand, the ethyl group continues to move away from the surface (see Figure 6). Atra interacts with Ca1 mostly in a bidentate way involving N1 and Cl. In the P2 snapshot, Ca1 – N1 = 2.548 Å and Ca1 – Cl = 3.132 Å. Let us notice that in the gas phase at B3LYP/6-31G*, those distances are slightly shorter, i.e.,
$Ca1 - N1 = 2.354 Å$ and $Ca1 - Cl = 2.878 Å$ [12,13]. The mean hydration of the pesticide is similar to M1, as well as the accessible $HS_{Ca1}$, whereas, $HS_{Ca2}$ increases from three to seven. The P2 snapshot shows the most frequent combination, namely $HS_{Ca1} = 1$ and $HS_{Ca2} = 6$, but $HS_{Ca1} = 2$ and $HS_{Ca2} = 5$ is almost equiprobable (see Table S1).

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Ce-surf, Ci-surf and Cl-surf normalized distance distributions for each of the singularities P1, M1, P2, M2 and P3 of the free energy profile and also for the $\zeta$ values 6.0, 7.0 and 8.0 Å.</figcaption></figure>

The biggest changes in inclination and hydration of Atra occur when moving from P2 to M2 corresponding to the largest free energy barrier of the profile (11.5 kcal·mol⁻¹). On average, six more water molecules enter the first HS of Atra (see Table 1). At the same time, $\bar{\theta}$ undergoes its biggest increase ever: +27.9° from 20.8 to 48.7°. Thus, Atra clearly quits its parallel orientation to the surface. Moreover, none of the alkyl groups conserve dispersion interactions with the surface. The angle increase and the loss of dispersion effects are obviously correlated and contribute to the high free energy barrier. The growth of $HS_{Atra}$ is not sufficient to compensate the main Atra-surf dispersion effects, which are lost. Only the Cl/surf dispersion remains, Cl-surf distance being always shorter than the reaction coordinate. The distribution presents two peaks centered on 2.4 Å and 3.6 Å. The chosen M2 snapshot shows 29 water molecules in $HS_{Atra}$, the clear remoteness of the alkyl groups as well as the Cl atom small distance to the surface with a $\theta$ value of 55.19°. The Ca1 hydration sphere in M2 includes zero to three molecules, with three being predominant (see Figure 7). At the same time, $HS_{Ca2}$ varies from one to four water molecules, with three being preferred (see Figure 8).
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Ca1-surf normalized distance distributions for each of the singularities P1, M1, P2, M2 and P3 of the free energy profile. The structures were sorted according to the number of water molecules in the first HS of Ca1 ($d_{Oe-Ca1} \le 3.0$ Å).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Ca2-surf normalized distance distributions for each of the singularities P1, M1, P2, M2 and P3 of the free energy profile. The structures were sorted according to the number of water molecules in the first HS of Ca2 ($d_{Oe-Ca2} \le 3.0$ Å).</figcaption></figure>
As a result of the small free energy difference between M2 and P3 (0.5 kcal·mol⁻¹), their descriptions are similar. When M2 presents several peaks with a normalized angle and distance distributions, in P3, only the most intense of those peaks is conserved (noted in bold for M2). For P3, the normalized distribution of θ is centered on 58° (M2 peaks around 31° and 59°). The greater the θ angle, the shorter the Cl-surf distance. Thus, in P3, the Cl-surf distance distribution is centered on 2.6 Å (Cl-surf M2 peaks around 2.4 Å and 3.6 Å). For those singularities, the shortest Cl-surf distances are found. Concerning the distributions of the alkyl group distance to the surface, Ce-surf and Ci-surf are centered on 6.8 Å and 6.6 Å, respectively (Ce-surf M2 peaks around 5.7 Å and 6.7 Å; Ci-surf M2 peaks around 5.5 Å and 6.5 Å). Concerning the cation hydration spheres, Ca1 and Ca2 are both surrounded by three water molecules. At the same time, HSₐtra varies from 24 to 35 for M2 and from 24 to 38 in P3, resulting in comparable mean values: 30.3 and 31.7, respectively. Finally, M2 resembles P3 more closely than its upstream minimum P2.

After P3, the free energy keeps on increasing without any clear singularity. To follow the evolution of the system when ζ is increased, three zones centered around 6.0, 7.0 and 8.0 Å were defined. For these zones, the majority of results are presented on Figures 3, 4 and 6 and on Tables 1 and 2. Cation-surf distance normalized distributions are presented in Figure S3 and the snapshots are depicted in Figure S4. The last interaction between Atra and the surface is through the Cl atom. This dispersion interaction is finally lost when Cl-surf neighbors 5.3 Å. In the 6.0 Å and 7.0 Å zones, there are structures that do or do not have this interaction, whereas in 8.0 Å zone, all the structures present long Cl-surf distances (>5.3 Å), no dispersion interaction between Atra and surf remains. The chosen snapshots (see Figure S4) shows this observation: in 6.0 Å and 7.0 Å, no water molecule lies in between the Cl atom of the pesticide and the Mont surface; for 8.0 Å, four water molecules enter the space between Cl and the surface. From P3, the Atra hydration sphere is growing but seems to reach a plateau after 7.0 Å with HSₐtra ≈ 36. In those cases, as Atra is going further to the surface, less water molecules are shared between Ca1 and Atra: 1.9 (over HSₐ1 = 3.0) around 7.0 Å and 1.1 (over HSₐ1 = 3.3) around 8.0 Å. Here, the maximum of HSₐtra includes 44 water molecules. As seen previously, a link can be made between θ angle distribution and Cl-surf distance distribution. In 6.0 Å, Cl-surf distribution presents two broad peaks: the first and most intense one centered on 3.2 Å corresponds to no water molecules between Cl and surf, simultaneously, the θ angle ranges from 25 to 69°; the second peak, with clearly less intensity, centered on 6.2 Å, gathers structures with more hydrated Cl, thus, the θ angle varies between -20 and +20°. Around 7.0 Å, the same tendency is observed but the peaks are closer for Cl-surf distance and θ angle distributions. When the Atra/surf interaction is lost for 8.0 Å, longer Cl-surf distances are accessible ([5.3; 7.8] Å). The corresponding angles range from 0° to 55° with a center value for this peak around 25°. A greater hydration of the Cl atom of the pesticide implies a smaller θ̄ value: 27.1° among the three last zones.

For all singularities and the three aforementioned zones, Ca1-Ca2 normalized distributions are presented in Figure S5, the corresponding mean Ca1-Ca2 distances are given in Table 1. As a reference, the Ca1-Ca2 distance in dry Mont is 10.341 Å [14,15]. In each case, the mean Ca1-Ca2 distance is shorter than or equal to the reference value. The greater the screening of the cation charge, the smaller the distance between Ca1 and Ca2. Screening is due to the presence of water molecules around both cations, but for Ca1, the pesticide also has a screening influence. The smallest distances are obtained for P1, M1 and P2 because of the Atra/Ca1 interaction. In M2 and P3, the interaction between Atra and Ca1 is lost. Despite larger HSₐ1 (2.1 in M2 and 3.0 in P3), the screening of Ca1 charge is weaker due to the loss of Atra/Ca1 interaction. Therefore, Ca1-Ca2 distances are larger than in P1, M1 and P2 but of the same magnitude as in dry Mont.
**3. Materials and Methods**

**3.1. Montmorillonite Model**

Montmorillonite is a 2:1 phyllosilicate containing one octahedral (O) sheet sandwiched between two tetrahedral (T) sheets. Oxygen (O²⁻) and hydroxyl (OH⁻) are located at the apices of the T and O polyhedrons. This dioctahedral smectite (T/O/T) has Si⁴⁺ cations in all sites of the T sheets and Al³⁺ in two thirds of the octahedral sites. Around 16.7% of the isomorphic substitutions of Al³⁺ by Mg²⁺ are distributed within the O layers. This distribution results in a net negative structural charge of the clay, which is balanced by adsorbed Ca²⁺ ions.

Unit cell parameters of pyrophyllite (a, b, c, α, β, γ) were found in good agreement with X-ray experimental data (R. Wardle et al. [16]) after being optimized at PBE-D2 level (see previous works [14,15,17]). Since the structure of pyrophyllite is similar to montmorillonite but without an isomorphic substitution nor a compensating cation in the interlayer spacing, it allowed us to build our montmorillonite model. The simulation cell consists of six unit cells (3a × 2b × c). The triclinic cell parameters are kept fixed during the Molecular Dynamics (MD) simulations at the following values: 15.500 Å × 17.931 Å × 25.000 Å and α = 91.030°, β = 100.374°, γ = 89.757°. The basal spacing *d*₀₀₁ is equal to 24.6 Å (*d*₀₀₁ = c sin β) resulting in an interlayer spacing δ ≈ 18 Å (see Figure 9).

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Ca²⁺–Montmorillonite (Mont) computational cell with the atrazine and 132 water molecules in the interlayer spacing. Mont contains four substitutions of Al³⁺ (light blue) by Mg²⁺ (yellow) in the octahedral (O) sheet and two Ca²⁺ cations (magenta) compensate the net negative charge carried by the clay layer and are situated above substitutions (ditrigonal ring of basal oxygen atoms, Ob). Ca1 is the Ca²⁺ interacting with Atra and Ca2 does not interact with Atra. The whole set of Ob atoms forms the Mont surface, noted surf. Si⁴⁺ ions in the tetrahedral (T) sheet are in dark blue. δ designates the interlayer spacing and *d*₀₀₁, the basal spacing. On the right, mechanism propositions for the desorption of Atra from Mont surface: (a) Atra alone is desorbed from Mont, Ca1 remains attached to Mont; (b) Atra–Ca1 complex is desorbed from Mont.</figcaption></figure>
## 3.2. Molecular Dynamics Simulations

MD simulations were performed using version 3.15.1 of CPMD software [18] within the Car-Parrinello framework [19] and a Density Functional Theory (DFT) potential using a generalized gradient approximation functional, PBE [20–23]. Grimme’s D2 corrections were selected to accurately model dispersion interactions [24–26]. The simulation cell was periodically replicated in 3D. Troullier-Martins norm-conserving pseudopotentials [27,28] with a plane-wave basis set describe valence-core interactions. A kinetic energy cutoff of 80 Ry for the plane-wave basis with a 4.0 a.u. time step was employed to integrate the equations of motion. Owing to the large computational cell, the Brillouin zone sampling was restricted to the Γ point using the Monkhorst-Pack method [29]. Considering the large size of the simulation cell, the default option for the Ewald term was used in the CPMD software. The fictitious mass of the electronic wavefunction was set to 500 a.u. The calculations were performed in a canonical (NVT) ensemble applying the Langevin thermostat procedure [30] at the simulation temperature of 350 K to avoid the over-structuration of water using PBE compared to experiment [31–33]. One hundred and thirty-two water molecules were present in the interlayer space in all the trajectories of this work.

## 3.3. Potential of Mean Force

Nineteen biased Car-Parrinello MD simulations were performed using one collective variable being our reaction coordinate noted ξ (see Table 3). The reaction coordinate ξ was defined as the distance, along the z-axis of the simulation cell, between the location of the center of mass of atrazine and the surface (surf) containing the basal oxygen atoms, Ob (see Figure 9). An umbrella sampling method [34–36] was used to observe the desorption of atrazine from the hydrated montmorillonite surface. Biased harmonic potentials centered on ξᵢ values were added to the system using PLUMED, an open-source community-developed library [37–39]. For a given ξᵢ, after an equilibration time of 4.0 ps, we performed NVT simulations lasting from 12 to 16 ps resulting in a cumulative simulation time of 209 ps.

<table><caption>**Table 3.** Parameters of the 19 biased MD trajectories<sup>a</sup>.</caption><thead><tr><th><i>i</i></th><th>ξ<sub>i</sub> (Å)</th><th><i>k</i><sub>i</sub> (kcal·mol<sup>-1</sup>·Å<sup>-2</sup>)</th><th><i>t</i><sub>i</sub> (ps)</th></tr></thead><tbody><tr><td>1</td><td>3.583</td><td>28.6</td><td>13</td></tr><tr><td>2</td><td>3.628</td><td>14.3</td><td>16</td></tr><tr><td>3</td><td>3.805</td><td>57.1</td><td>4</td></tr><tr><td>4</td><td>3.805</td><td>28.6</td><td>9</td></tr><tr><td>5</td><td>4.157</td><td>14.3</td><td>12</td></tr><tr><td>6</td><td>4.486</td><td>14.3</td><td>14</td></tr><tr><td>7</td><td>5.216</td><td>14.3</td><td>12</td></tr><tr><td>8</td><td>5.421</td><td>57.1</td><td>5</td></tr><tr><td>9</td><td>5.421</td><td>28.6</td><td>7</td></tr><tr><td>10</td><td>5.615</td><td>28.6</td><td>13</td></tr><tr><td>11</td><td>5.745</td><td>14.3</td><td>14</td></tr><tr><td>12</td><td>6.274</td><td>14.3</td><td>12</td></tr><tr><td>13</td><td>6.803</td><td>14.3</td><td>14</td></tr><tr><td>14</td><td>7.332</td><td>14.3</td><td>14</td></tr><tr><td>15</td><td>7.832</td><td>57.1</td><td>4</td></tr><tr><td>16</td><td>7.832</td><td>28.6</td><td>10</td></tr><tr><td>17</td><td>7.861</td><td>14.3</td><td>12</td></tr><tr><td>18</td><td>8.391</td><td>14.3</td><td>12</td></tr><tr><td>19</td><td>8.920</td><td>14.3</td><td>12</td></tr></tbody></table>

Spanning the range from 3.583 Å to 8.920 Å for the reaction coordinate ξ allowed to move atrazine from a location close to the surface to the middle of the interlayer spacing. The sampling interval changed according to the free energy profile singularities

<footnote><sup>a</sup> *i* is the window index, ξᵢ the value of the reaction coordinate around which the harmonic potential oscillated, *k*ᵢ is the associated harmonic constant, and *t*ᵢ the duration of the simulation (ps).</footnote>
(from 0.029 Å minimum to 0.529 Å maximum) in order to obtain satisfying overlapping distributions of $\zeta$ (see Figure S1).

The Potential of Mean Force (PMF) was derived to quantify the free energy $F$ of the atrazine desorption from hydrated montmorillonite surface. The changes in $F$ were due to the changes in the system configurations and were described by the PMF [40–43]. In an NVT statistical ensemble, the free energy $F$ of a system in thermodynamic equilibrium could be written as:

$$F = -k_B T \ln Z \qquad (1)$$

$Z$ being the canonical partition function, $T$, the temperature (K) and $k_B$, the Boltzmann constant ($J \cdot K^{-1}$).

For each $i$ window, the dynamic simulation is biased introducing into the Hamiltonian of the system, a harmonic potential $V$ centered on the $\zeta_i$ value:

$$V = \frac{1}{2} k_i (\zeta - \zeta_i)^2 \qquad (2)$$

$k_i$ is the harmonic constant (see Table 3) for window $i$. The values of the $k_i$ were selected to obtain the best sampling of the potential energy surface (PES).

For the hydrated interacting system under study (Atra;Mont), a free energy curve was obtained for each window $i$ noted $W(\zeta)$ being the PMF and expressed by the following equation:

$$W(\zeta) = -k_B T \ln(g(\zeta)) \qquad (3)$$

$g(\zeta)$ was the distribution function of distances for the pair (Atra; Mont).

The entire free energy curve was rebuilt from the $\zeta$ distributions of each biased MD trajectory via the WHAM self consistent method (Weighted Histogram Approximation Method) [41,44,45]. The WHAM parameters were set to 1000 for the number of bins and to $10^{-6}$ kcal·mol⁻¹ for the precision on final free energy. Thus, the number of points in the free energy curve is 1000 and the WHAM iterative resolution process is considered to be converged when no free energy value for any simulation window changes by more than $10^{-6}$ kcal·mol⁻¹ on consecutive iterations.

## 4. Conclusions

In the beginning of the desorption process, Atra interacts with the surface through Ca1. The first two minima, P1 and P2, are free energetically quasi degenerated. In P1, Atra bonds to Ca1 via the monodentate site N1, whereas in P2, both N1 and Cl are involved in the interaction resulting in a bidentate binding site. In P1, Atra is parallel to the Mont surface maximizing the Atra/surf dispersion interactions. Moving to P2, Atra changes its orientation from parallel to inclined. The M1 barrier connecting P1 to P2 is low in free energy (≈0.5 kcal·mol⁻¹) and corresponds to the remoteness of the ethyl group with respect to the surface and thus, reducing Atra/surf dispersion effects. Between P2 and M2, the barrier is clearly higher (≈11.5 kcal·mol⁻¹). This is due to the loss of isopropyl/surf interaction, which is mainly dispersive, for Atra only the Cl/surf dispersion remains. M2 resembles more P3 than its upstream minimum P2. From P3, when $\zeta$ keeps on increasing, the free energy also increases without any clear singularity. Finally, around $\zeta = 8.0$ Å, the last Atra/surf dispersion interaction via the Cl atom vanishes and Atra is completely desorbed. The mechanism evidenced in the desorption process is the one presented in Figure 9a: Atra desorbed alone and Ca1 remaining on the Mont surface, which is similar to the previously studied fenhexamid (Fen) pesticide desorption from the same clay surface [11]. However, the Atra desorption profile differs from the Fen desorption profile. First, the highest barrier to desorb Fen from Mont is ≈6.5 kcal·mol⁻¹, which is about half the barrier to desorb Atra. Moreover in the case of Fen, the desorption is driven by both cation hydration, whereas here, the process is led by the dispersion effects of the Atra alkyl groups. As the mechanism is driven by atrazine alkyl groups and their sizes, we can deduce the desorption behavior of the various atrazine metabolites. Therefore, we project an easier desorption for
deisopropylatrazine, than atrazine or deethylatrazine. It would be interesting to confirm this hypothesis by simulating the desorption processes of Atra degradation products.

**Supplementary Materials:** The supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms25031604/s1.

**Author Contributions:** Conceptualization, F.B. and S.H.; methodology, F.B. and S.H.; calculation, Q.D.; formal analysis, F.B. and S.H.; investigation, Q.D.; resources, F.B. and S.H.; writing—original draft preparation, Q.D.; writing—review and editing, F.B. and S.H.; funding acquisition, F.B. and S.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was co-funded by Région Occitanie, France and Ecole d'Ingénieurs de Purpan, grant number [ADSORCO-21013061].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Additional data are already joined in supporting information.

**Acknowledgments:** This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2022-[P1222] and to the HPC resources of TGCC under the allocation 2022-[cpq7002] made by GENCI.

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

## References

1. Vieira, Y.; Silveira, J.P.; Dotto, G.L.; Knani, S.; Vieillard, J.; Georgin, J.; Franco, D.S.P.; Lima, E.C. Mechanistic insights and steric interpretations through statistical physics modelling and density functional theory calculations for the adsorption of the pesticides atrazine and diuron by Hovenia dulcis biochar. *J. Mol. Liq.* **2022**, *367*, 120418. [CrossRef]

2. Cheng, Y.; Wang, B.; Shen, J.; Yan, P.; Kang, J.; Wang, W.; Bi, L.; Zhu, X.; Li, Y.; Wang, S.; et al. Preparation of novel N-doped biochar and its high adsorption capacity for atrazine based on Pi-Pi electron donor-acceptor interaction. *J. Hazard. Mater.* **2022**, *432*, 128757. [CrossRef] [PubMed]

3. Vieira, Y.; Schnorr, C.; Piazzi, A.C.; Netto, M.S.; Piccini, W.M.; Franco, D.S.P.; Mallmann, E.S.; Georgin, J.; Silva, L.F.O.; Dotto, G.L. An advanced combination of density functional theory simulations and statistical physics modeling in the unveiling and prediction of adsorption mechanisms of 2,4-D pesticide to activated carbon. *J. Mol. Liq.* **2022**, *361*, 119639. [CrossRef]

4. Jablonowski, N.D.; Schäffer, A.; Burauel, P. Still present after all these years: Persistence plus potential toxicity raise questions about the use of atrazine. *Environ. Sci. Pollut. Res.* **2011**, *18*, 328–331. [CrossRef] [PubMed]

5. Council Regulation (EEC) No. 793/93 on the Evaluation and Control of the Risks of Existing Substances. *Official Journal L* 84. 5 April 1993; pp. 1–75. Available online: https://www.fao.org/faolex/results/details/en/c/LEX-FAOC037952 (accessed on 21 January 2024).

6. Rostami, S.; Jafari, S.; Moeini, Z.; Jaskulak, M.; Keshtgar, L.; Badeenezhad, A.; Azhdarpoor, A.; Rostami, M.; Zorena, K.; Dehghani, M. Current methods and technologies for degradation of atrazine in contaminated soil and water: A review. *Environ. Technol. Innov.* **2021**, *24*, 102019. [CrossRef]

7. Abate, G.; Masini, J.C. Adsorption of Atrazine, Hydroxyatrazine, Deethylatrazine, and Deisopropylatrazine onto Fe(III) Polyhydroxy Cations Intercalated Vermiculite and Montmorillonite. *J. Agric. Food Chem.* **2005**, *53*, 1612–1619. [CrossRef]

8. Han, Y.; Gu, L.; Zhang, M.; Li, Z.; Yang, W.; Tang, X.; Xie, G. Computer-aided design of molecularly imprinted polymers for recognition of atrazine. *Comput. Theor. Chem.* **2017**, *1121*, 29–34. [CrossRef]

9. França, V.L.B.; Amaral, J.L.; Martins, Y.A.; Caetano, E.W.S.; Brunaldi, K.; Freire, V.N. Characterization of the binding interaction between atrazine and human serum albumin: Fluorescence spectroscopy, molecular dynamics and quantum biochemistry. *Chem.-Biol. Interact.* **2022**, *366*, 110130. [CrossRef]

10. Moreno-Rodriguez, D.; Jankovic, L.; Scholtzova, E.; Tunega, D. Stability of Atrazine–Smectite Intercalates: Density Functional Theory and Experimental Study. *Minerals* **2021**, *11*, 554. [CrossRef]

11. Belzunces, B.; Hoyau, S.; Cuny, J.; Bessac, F. Desorption Mechanism of a Pesticide from a Hydrated Calcium Montmorillonite Unraveled by Molecular Dynamics Simulations. *J. Phys. Chem. C* **2023**, *127*, 12953–12966. [CrossRef]

12. Bessac, F.; Hoyau, S. Pesticide interaction with environmentally important cations: A theoretical study of atrazine. *Comput. Theor. Chem.* **2011**, *966*, 284–298. [CrossRef]

13. Bessac, F.; Hoyau, S. Pesticide interaction with environmentally important cations: A theoretical study of atrazine in interaction with two Ca²⁺ cations. *Comput. Theor. Chem.* **2013**, *1022*, 6–13. [CrossRef]

14. Belzunces, B.; Hoyau, S.; Benoit, M.; Tarrat, N.; Bessac, F. Theoretical study of the atrazine pesticide interaction with pyrophyllite and Ca²⁺-montmorillonite clay surfaces. *J. Comput. Chem.* **2016**, *38*, 133–143. [CrossRef] [PubMed]
15. Belzunces, B.; Hoyau, S.; Bessac, F. Interaction of Metamitron and Fenhexamid with Ca²⁺-Montmorillonite Clay Surfaces: A Density Functional Theory Molecular Dynamics Study. *J. Comput. Chem.* **2019**, *40*, 1449–1462. [CrossRef] [PubMed]

16. Wardle, R.; Brindley, G.W. The crystal structures of pyrophyllite, 1Tc, and of its dehydroxylate. *Am. Mineral.* **1972**, *57*, 732–750.

17. Belzunces, B.; Hoyau, S.; Cuny, J.; Bessac, F. Pesticide interaction with environmentally important cations: A molecular dynamics and DFT study of metamitron and fenhexamid. *Comput. Theor. Chem.* **2017**, *1117*, 220–234. [CrossRef]

18. CPMD. Copyright IBM Corp 1990–2015, Copyright MPI für Festkörperforschung Stuttgart 1997–2001. Available online: https://github.com/CPMD-code (accessed on 21 January 2024).

19. Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. *Phys. Rev. Lett.* **1985**, *55*, 2471–2474. [CrossRef]

20. Kohn, W.; Sham, L.J. Self-Consistent Equations Including Exchange and Correlation Effects. *Phys. Rev.* **1965**, *140*, A1133–A1138. [CrossRef]

21. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868; Erratum in *Phys. Rev. Lett.* **1997**, *78*, 1396. [CrossRef]

22. Voora, V.K.; Al-Saidi, W.A.; Jordan, K.D. Density Functional Theory Study of Pyrophyllite and M-Montmorillonites (M = Li, Na, K, Mg, and Ca): Role of Dispersion Interactions. *J. Phys. Chem. A* **2011**, *115*, 9695–9703. [CrossRef]

23. Tunega, D.; Bučko, T.; Zaoui, A. Assessment of ten DFT methods in predicting structures of sheet silicates: Importance of dispersion corrections. *J. Chem. Phys.* **2012**, *137*, 114105. [CrossRef]

24. Grimme, S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. *J. Comput. Chem.* **2004**, *25*, 1463–1473. [CrossRef]

25. Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. *J. Comput. Chem.* **2006**, *27*, 1787–1799. [CrossRef]

26. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. *J. Chem. Phys.* **2010**, *132*, 154104. [CrossRef]

27. Hamann, D.R.; Schlüter, M.; Chiang, C. Norm-Conserving Pseudopotentials. *Phys. Rev. Lett.* **1979**, *43*, 1494–1497. [CrossRef]

28. Troullier, N.; Martins, J.L. Efficient pseudopotentials for plane-wave calculations. *Phys. Rev. B* **1991**, *43*, 1993–2006. [CrossRef]

29. Monkhorst, H.J.; Pack, J.D. Special points for Brillouin-zone integrations. *Phys. Rev. B* **1976**, *13*, 5188–5192. [CrossRef]

30. Ceriotti, M.; Bussi, G.; Parrinello, M. Langevin Equation with Colored Noise for Constant-Temperature Molecular Dynamics Simulations. *Phys. Rev. Lett.* **2009**, *102*, 020601. [CrossRef] [PubMed]

31. Soper, A.K. The radial distribution functions of water and ice from 220 to 673 K and at pressures up to 400 MPa. *Chem. Phys.* **2000**, *258*, 121–137. [CrossRef]

32. Skinner, L.B.; Huang, C.; Schlesinger, D.; Pettersson, L.G.M.; Nilsson, A.; Benmore, C.J. Benchmark oxygen-oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range. *J. Chem. Phys.* **2013**, *138*, 074506. [CrossRef] [PubMed]

33. Bankura, A.; Karmakar, A.; Carnevale, V.; Chandra, A.; Klein, M.L. Structure, Dynamics, and Spectral Diffusion of Water from First-Principles Molecular Dynamics. *J. Phys. Chem. C* **2014**, *118*, 29401–29411. [CrossRef]

34. Kirkwood, J.G. Statistical Mechanics of Fluid Mixtures. *J. Chem. Phys.* **1935**, *3*, 300–313. [CrossRef]

35. Torrie, G.M.; Valleau, J.P. Monte Carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical Lennard-Jones fluid. *Chem. Phys. Lett.* **1974**, *28*, 578–581. [CrossRef]

36. Torrie, G.M.; Valleau, J.P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. *J. Comput. Phys.* **1977**, *23*, 187–199. [CrossRef]

37. The PLUMED Consortium. Promoting transparency and reproducibility in enhanced molecular simulations. *Nat. Methods* **2019**, *16*, 670. [CrossRef] [PubMed]

38. Tribello, G.A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED2: New feathers for an old bird. *Comp. Phys. Comm.* **2014**, *185*, 604. [CrossRef]

39. Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R.A.; et al. PLUMED: A portable plugin for free energy calculations with molecular dynamics. *Comp. Phys. Comm.* **2009**, *180*, 1961. [CrossRef]

40. Kollman, P. Free energy calculations: Applications to chemical and biochemical phenomena. *Chem. Rev.* **1993**, *93*, 2395–2417. [CrossRef]

41. Kumar, S.; Rosenberg, J.M.; Bouzida, D.; Swendsen, R.H.; Kollman, P.A. Multidimensional free-energy calculations using the weighted histogram analysis method. *J. Comput. Chem.* **1995**, *16*, 1339–1350. [CrossRef]

42. Roux, B. The calculation of the potential of mean force using computer simulations. *Comput. Phys. Commun.* **1995**, *91*, 275–282. [CrossRef]

43. Trzesniak, D.; Kunz, A.P.E.; van Gunsteren, W.F. A Comparison of Methods to Compute the Potential of Mean Force. *ChemPhysChem* **2007**, *8*, 162–169. [CrossRef] [PubMed]
44. Kumar, S.; Rosenberg, J.M.; Bouzida, D.; Swendsen, R.H.; Kollman, P.A. The weighted histogram analysis method for free-energy calculations on biomolecules. *I*. The method. *J. Comput. Chem.* **1992**, *13*, 1011–1021. [CrossRef]

45. Grossfield, A. WHAM: *The Weighted Histogram Analysis Method*, Version 2.0.9.1. Available online: http://membrane.urmc.rochester.edu/wordpress/?page_id=126 (accessed on 21 January 2024).

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