minerals Article Adsorption of Pharmaceuticals onto Smectite Clay Minerals: A Combined Experimental and Theoretical Study Gwenaëlle Corbin 1 , Emmanuelle Vulliet 2 , Bruno Lanson 3 , Albert Rimola 4, * and Pierre Mignon 1, * 1 Institut Lumière Matière, UMR5306 Université Claude Bernard Lyon 1-CNRS, Université de Lyon, 69622 Villeurbanne CEDEX, France; gw.corbin@hotmail.fr 2 Univ Lyon, CNRS, Université Claude Bernard Lyon 1, Institut des Sciences Analytiques, UMR5280, 5 rue de la Doua, 69100 Villeurbanne, France; emmanuelle.vulliet@isa-lyon.fr 3 Univ. Grenoble Alpes, Univ. Savoie Mont-Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, F-38000 Grenoble, France; bruno.lanson@univ-grenoble-alpes.fr 4 Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalonia, Spain * Correspondence: albert.rimola@uab.cat (A.R.); pierre.mignon@univ-lyon1.fr (P.M.); Tel.: +33-0472448314 (P.M.) Abstract: The adsorption of two pharmaceuticals, carbamazepine and paracetamol, onto the expand- able clay mineral saponite has been studied through the combination of kinetic experiments, X-ray diffraction, and theoretical modeling. Kinetic experiments indicate low adsorption for carbamazepine and paracetamol on expandable smectite clay. Accordingly, X-ray diffraction experiments show that neither compound enters smectite interlayer space. Molecular dynamics simulations were carried out to understand the interactions between the two pharmaceuticals and the saponite basal surface in the presence of Na+ cations. Calculations reveal that paracetamol almost does not coordinate solution cations, whereas a rather low coordination to cation is observed for carbamazepine. As a result, the adsorption onto the clay surface results mainly from van der Waals interactions for both pharmaceuticals. Carbamazepine does adsorb the surface via two configurations, one involving   cation coordination, which corresponds to a rather stable adsorption compared to paracetamol. This Citation: Corbin, G.; Vulliet, E.; is confirmed by structural analyses completed with desorption free energy profile. Lanson, B.; Rimola, A.; Mignon, P. Adsorption of Pharmaceuticals onto Keywords: organic contaminants; paracetamol; carbamazepine; smectite; interlayer space; molecular Smectite Clay Minerals: A Combined dynamics; kinetic experiments; X-ray diffraction Experimental and Theoretical Study. Minerals 2021, 11, 62. https://doi. org/10.3390/min11010062 1. Introduction Received: 30 November 2020 For about thirty years, the concern about the amount of organic contaminants that Accepted: 7 January 2021 can be found in various environmental compartments has grown [1–3], as massively used Published: 11 January 2021 pharmaceutical compounds (PCs) have been continuously introduced in the environment over time. In addition to unused medication and pharmaceutical manufacturing plants and Publisher’s Note: MDPI stays neu- tral with regard to jurisdictional clai- hospitals discharge, a large fraction of consumed human or veterinary active compounds ms in published maps and institutio- are expelled from the body after a short period of time and end up in municipal sewage nal affiliations. systems and wastewater treatment plants [4–9]. Land applications of treated sewage sludge contribute to nutrient recycling and soil fertilization but also lead to increasing PC concentrations in soils such as agricultural fields as PCs do not always degrade during wastewater treatments [10–12], nor in soils [13,14]. Retention of persistent PCs in sludge Copyright: © 2021 by the authors. Li- or in soil minerals potentially prevents their release in groundwater where they can re- censee MDPI, Basel, Switzerland. enter the biosphere and drinking water supplies [15,16]. As such, the adsorption of This article is an open access article PCs by minerals is an important issue either to predict the PCs’ soil retention properties distributed under the terms and con- or to develop minerals which, when added to sludge, may retain PCs, thus enhancing ditions of the Creative Commons At- biodegradation before land spreading. tribution (CC BY) license (https:// Several studies investigated the degradation and mobility of persistent PCs (e.g., creativecommons.org/licenses/by/ carbamazepine, diclofenac, or ibuprofen) in soils as a function of their mineralogy and/or 4.0/). Minerals 2021, 11, 62. https://doi.org/10.3390/min11010062 https://www.mdpi.com/journal/minerals Minerals 2021, 11, 62 2 of 12 organic matter composition. PC degradation is directly related to the presence of organic matter, which also reduces their mobility [13,14,17]. The amount of clays was also observed to have an impact, although less important than that of organic matter, on soil retention properties by decreasing PC mobility. Retention studies were also performed as a function of soil composition (clayey/sandy loam) for various pharmaceuticals and steroid hormones, which also confirmed the role of clay materials in PC retention [18,19]. Clays minerals are abundant in nature and especially relevant for their adsorptive properties owing to their layered structure and ionic exchange properties. They have thus been widely used as environmental decontamination and low-cost substitution agents and as materials for the removal of heavy metals from wastewaters [20,21]. Particular interest is devoted also to their ability to adsorb organic compounds, especially pharmaceuti- cals [22]. Recent studies investigated in detail the adsorption of PCs on sodium-exchanged montmorillonite under various pH conditions [23]. The three investigated PCs (codeine, diazepam, and oxazepam) were found to be intercalated into the interlayer space via cation exchange at acidic pH, while at neutral pH, only codeine was protonated and present in clay interlayers. Among a large pool of PCs, it was also found that anionic PCs may adsorb to clay edges, while neutral ones may interact through van der Waals interactions [24]. Theoretical calculations are useful tools to understand further the adsorption of or- ganic compounds onto clay surfaces, and they have already shown their ability to provide a precise image of the molecular interactions at play in the adsorption mechanism [25–28]. For example, modeling the adsorption of two glyphosate ionic forms on kaolinite sur- face [29] showed a larger adsorption free energy for the anionic form. Another study coupling X-ray diffraction (XRD), infrared and nuclear magnetic resonance spectroscopies, and molecular dynamics (MD) simulations provided a detailed description of oxytetracy- cline interactions with sodium-montmorillonite [30]. MD simulations indicated that clay charge localization controlled the binding conformation of intercalated oxytetracycline, allowing the formation of multiple interactions consistent with the spectroscopic data. In the present study, kinetic experiments combined to XRD measurements and MD simulations are used to describe the adsorption of two PCs (carbamazepine and parac- etamol) onto a smectite clay (saponite) and to understand the mechanisms driving these sorptive interactions. We will focus on these two PCs because different behaviors in soil were observed in previous studies [18]. 2. Materials and Methods 2.1. Materials Saponite ([Na0.8 ]inter [Mg6.0 ]octa[Si7.2 Al0.8 ]tetra O20 (OH)4 ; hereafter referred to as Sap) was synthesized hydrothermally from a gel precursor of adequate stoichiometry [31]. Syn- thesis was performed in an externally heated Morey-type pressure vessel with an internal silver tubing [32,33]. Synthesis conditions were a temperature of 400 ◦ C, a water pressure of 1 kbar, and a duration of 4 weeks, as described elsewhere [31]. After synthesis, the Sap was Na-saturated by contact with a 1 mol L–1 aqueous NaCl solution with mechanical shaking for 24 h to ensure a complete exchange of interlayer cations before separation of the solid fraction by centrifugation. Excess NaCl was removed by washing the solid three times with deionized water (Siemens UltraClear, Metagenics, Inc., Aliso Viejo, CA, USA 18.2 MΩ cm–1 ) and separation of the solid fraction. Sap cation exchange capacity was previously measured at 105 meq/100 g [34]. 2.2. Kinetic Experiments Low PC concentrations have been considered in order to reflect values observed in the environment. The presence of PCs in surface or ground water has been widely studied, in which carbamazepine and, to a lesser extent, paracetamol rank among the most often found molecules [35–37]. In the extensive monitoring they carried out across Europe, Loos et al. reported carbamazepine in 42% of analyzed groundwaters (maximum concentration 0.39 µg/L) and in 95% of surface waters (maximum concentration 11.561 µg/L, average Minerals 2021, 11, 62 3 of 12 0.075 µg/L). Paracetamol is found at levels ranging from a few ng/L to 1.5 µg/L in water resources [37]. To follow the adsorption kinetics of the two PCs onto Sap, 200 mg of clay was introduced into tubes containing paracetamol or carbamazepine solutions (V = 1 mL, C = 2 µg/L). This mixture was then stirred in the dark at room temperature. Tubes were removed at regular times, centrifuged, and the supernatant was analyzed by liquid chromatography coupled to tandem mass spectrometry (LC-MS/MS). The amount of PC adsorbed onto the clay was determined indirectly as the difference between the quantity initially introduced into the tube and the quantity remaining at the time of sampling. We would like to stress at that such low concentrations (consistent with those found in natural environments), interlayer sorption is expected to be unlikely. Using higher concen- trations would likely involve different surface sites and thus different, environmentally non-relevant sorption processes. In addition, one may note that even at such low concen- trations, sorption is not complete, supporting the low affinity of the two PCs investigated for smectite basal surfaces (and interlayers). Analyses were performed using LC-MS/MS with a chromatograph system Agilent 1290 series liquid (Agilent Technologies, Avondale, AZ, USA) coupled to a 3200 QTrap (AB Sciex, CA, USA) triple quadrupole mass spectrometer equipped with an electro- spray ion source (ESI Turbo V) operated in positive mode. A Zorbax Eclipse C18 column (50 mm × 2.1 mm i.d., 1.8 µm) from Agilent was used for the chromatographic separation, set at 60 ◦ C. The mobile phase was composed of Milli-Q quality water with 0.01% acetic acid and methanol. The flow rate was 0.6 mL/min, and the sample volume injected was 50 µL. The quantification of paracetamol and carbamazepine was achieved in multiple re- action monitoring mode (MRM). The two MRM transitions 151.8 → 110.2 and 151.8 → 64.9 were used for paracetamol; 236.9 → 193.9.2 and 236.9 → 192 were used for carbamazepine. Nitrogen was used as the drying and nebulizing gas. All data were collected and processed using Analyst 1.5.2 software. 2.3. X-ray Diffraction For all samples, oriented slides were prepared by drying a clay suspension on glass slides at room temperature following contact with the organic solution. X-ray diffraction (XRD) patterns were then recorded using a Bruker D8 diffractometer equipped with an MHG Messtechnik humidity controller coupled to an Anton Paar CHC+ chamber. Intensi- ties were measured with a SolXE Si(Li) solid-state detector (Baltic Scientific Instruments) for 4 s per 0.04◦ 2θ step over the 2−50◦ 2θ Cu Kα angular range. The divergence slit, the two Soller slits, the antiscatter, and resolution slits were 0.3◦ , 2.3◦ , 0.3◦ , and 0.1◦ , respectively. Samples were kept at 23 ◦ C in the CHC+ chamber during data collection. A constant flow of mixed dry/saturated air maintained a 40% relative humidity (RH) during data collection after an initial equilibration at ∼97% RH. RH was continuously monitored with a hygrometer (uncertainty of ∼2% RH) located close to the sample. 2.4. Computational Details The saponite layer consisted of two tetrahedral (T) SiO4 sheets sandwiching an octahe- dral (O) MgO6 (OH)2 one. As saponite is trioctahedral, all octahedral sites are occupied, and only tetrahedral Si4+ atoms are substituted with Al3+ , leading to a negative layer charge compensated for by exchangeable counter ions present in the interlayer space. The general formula for the saponite studied here is (Si7.2Al0.8 )(Mg6 )O20 (OH)4Na0.8 , which thus has a 0.8e charge per unit cell. An orthogonal model was considered with a 6 × 4 supercell with dimensions of a = 32.22 Å and b = 37.12 Å within the layer plane and c = 60 Å. The con- sidered cell included 19 tetrahedral substitutions and, thus, 19 Na+ counter ions. Atomic substitutions were randomly located, respecting, however, a separation rule between substitutions sites of at least one tetrahedral center between two substitutions. In the simulation box, single molecules of the PCs were introduced. Both carba- mazepine and paracetamol were considered in their neutral forms because of their high Minerals 2021, 11, 62 4 of 12 pKa. Geometries and charge parameters were refined through DFT calculations at the B3LYP/6–31G* level. Molecules were optimized in the gas phase, in which the restrained electrostatic potential (RESP) partial atomic charges reproducing the electrostatic poten- tial of the molecular boundaries were used for initial guess adsorptions [38]. It is worth bearing in mind that these simulated systems are model systems, exclusively focused on investigating a particular problem: the interaction of PC with clays. Hence, we use one PC molecule interacting with a smectite clay surface model, which includes most of its surface features. The use of more molecules in the model would not add more information to the studied problem. The aim of the simulations is to understand the sorption mechanism and to explain it through the analysis (strength and type) of the molecular interactions taking part in the sorption. The ClayFF force field was used for the clay model [39], whereas water and ions were modeled using the extended simple point charge (SPC/E) and Dang models, respec- tively [40,41]. The SPC/E model is known to accurately reproduce water’s structure [42]. The general Amber force field [43] was used to generate the PCs’ molecular parameters and Lorentz–Berthelot rules were applied for inter-force fields’ parameters. This approach has already been successfully used in various theoretical works [28,44–46]. The Amber force field has been extensively used for organic molecules’ modeling and reproduces well their structure and dynamics. Amber and ClayFF are rather simple to combine as they both use harmonic potential for bond terms, although ClayFF does not use angle and dihedral terms. Periodic boundary conditions were applied in all three spatial directions. The large-scale atomistic/molecular massively parallel simulator (LAMMPS) [47] was used for all simulations. A unique clay layer was considered to model the clay basal surface, divided with half of it at the bottom and at the top of the box along the vertical parameter. The initial separation between the lower and upper surfaces was set to 60 Å, with PC molecules being placed at the box’s center. Water molecules were included to fill the empty space using a regular grid (with an interval of 3.3 Å) and including the amount of Na+ needed to ensure charge compensation. A cutoff of 12 Å was applied for the electrostatic and Lennard–Jones interactions. Coulombic interactions were computed using the Ewald summation and the particle-particle-particle-mesh method at a 0.0001 accuracy. All simulations were carried out at 300 K and 1 atm. The oxygen–hydrogen bonds of water molecules were kept fixed by applying the shake algorithm [48]. The equilibration period consisted of warming up the system, starting from water/ions moiety, then the clay particles, and finally, PC molecules in the canonical NVT ensemble. The system was then allowed to relax along all directions through the isothermal–isobaric NPT ensemble to model the hydrated interlayer space in real conditions and density. The vertical separation between basal surfaces was then kept fixed during the whole production time. The system was then equilibrated through a two-step procedure, increasing the time step from 0.1 to 0.5 fs. The whole equilibration run included a total of 800,000 time steps. During the whole simulation and production, octahedral Mg atoms were kept fixed. The production time during which data were collected and analyzed was 100 ns. The desorption free energy profiles were obtained through umbrella sampling [49] via the plumed module [50,51] patched to LAMMPS. These simulations were performed using a series of productions where the molecule–surface distance (along z) is varied, from an adsorbed configuration to a desorbed one. The vertical distance was computed between the center of benzene cycles of the molecules to the average surface oxygen atoms for the clay surface. PC molecules were thus allowed to rotate but kept parallel to the surface. Each simulation step of the umbrella sampling consisted of 4 million steps (∆t = 0.5 fs). Molecule–surface distances were varied from 2.5 to 8.5 Å with a 0.1-Å displacement step (61 simulations were thus performed for each molecule). The potential of mean force (PMF) was integrated over all simulations. Minerals 2021, 11, 62 5 of 12 Minerals 2021, 11, x FOR PEER REVIEW 5 of 13 3. Results 3. Results 3.1. Kinetics 3.1. Kinetics As shown As shown in in Figure Figure 1, 1, the the adsorption adsorption ofof both both carbamazepine carbamazepine andand paracetamol paracetamol ontoonto clay is clay is very veryfast fastover overthe thefirst hours first of of hours contact time, contact andand time, maximal retention maximal is reached retention after is reached the initial after first hours. the initial After reaching first hours. the saturation After reaching value in value the saturation about in 50 about hours 50 of contact time, h of contact the quantity time, adsorbed the quantity no longer adsorbed varies significantly. no longer The small varies significantly. Theamounts adsorbed small amounts at equi- adsorbed at equilibrium (35% and 20% for carbamazepine and paracetamol, respectively) arenote- librium (35% and 20% for carbamazepine and paracetamol, respectively) are also also worthy. noteworthy. Figure 1. Sorption Figureof1.carbamazepine (a) and paracetamol Sorption of carbamazepine (b) onto saponite (a) and paracetamol claysaponite (b) onto at 200 g/L. clay at 200 g/L. 3.2. X-ray Diffraction (XRD) For the two series of samples, the XRD patterns exhibit a series of basal reflections reflections atat ~14.9, 7.7, 5.07, 3.85, and 3.05 Å (5.9, 11.45, 17.5, 23.4, and 29.2 ◦ 2θ Cu Kα, respectively), ~14.9, 7.7, 5.07, 3.85, and 3.05 Å (5.9, 11.45, 17.5, 23.4, and 29.2 °2θ Cu Kα, respectively), typical typical of of bi-hydrated bi-hydratedsmectites smectites(see (seeFigure Figure22) )[52,53]. [52,53].The Theslight slightdeparture departurefrom from rationality rational- indicates the coexistence ity indicates the coexistenceof different hydration of different statesstates hydration and theandinterstratification of a minor the interstratification of a amount of mono-hydrated minor amount layers interstratified of mono-hydrated in mainly in layers interstratified bi-hydrated smectite crystals. mainly bi-hydrated In smectite addition, a shoulder is visible at ~12.8 Å (~6.9 ◦ 2θ Cu Kα) on the high-angle side of the first crystals. In addition, a shoulder is visible at ~12.8 Å (~6.9 °2θ Cu Kα) on the high-angle basal side ofreflection, the first indicative of minor basal reflection, domainsof indicative in minor which domains mono-hydrated in whichlayers prevail. The mono-hydrated dual hydration state exhibited by the samples is likely related to the presence of bothtoCa 2+ layers prevail. The dual hydration state exhibited by the samples is likely related the (from solution during the sorption experiments) + and Na experiments) (from the initial clay presence of both Ca2+ (from solution during the sorption and Nasaturation) + (from the cations in saponite interlayers. initial clay saturation) cations in saponite interlayers. In any case, the XRD patterns for all samples collected during the sorption isotherms are essentially similar, showing no significant evolution with increasing contact time with carbamazepine or paracetamol. This remarkable stability of the layer-to-layer distance indicates that neither paracetamol nor carbamazepine molecules enter the Sap interlayer space, consistent with the low sorption capacity determined from kinetic experiments. Rather, Ca2+ and Na+ cations persist in these interlayers to compensate for the layer charge deficit, thus keeping constant the overall hydration state of Sap at a given RH. Minerals 2021, 11, 62 6 of 12 Minerals 2021, 11, x FOR PEER REVIEW 6 of 13 a) blank exp. Carbamazepine - 1h Carbamazepine - 2h Carbamazepine - 3h Carbamazepine - 5h Carbamazepine - 7h Carbamazepine - 8h Carbamazepine - 24h 2 4 6 8 10 11 14 17 20 23 26 29 32 35 38 41 44 47 50 b) blank exp. Paracetamol - 1h Paracetamol - 2h Paracetamol - 3h Paracetamol - 5h Paracetamol - 7h Paracetamol - 8h Paracetamol - 24h 2 4 6 8 10 11 14 17 20 23 26 29 32 35 38 41 44 47 50 Position (°2θ CuKα) Figure 2. X-ray Figure 2. X-ray diffraction diffraction patterns patterns of of saponite saponite (Sap) (Sap) contacted contacted with with (a) carbamazepine and (a) carbamazepine and (b) (b) paracetamol paracetamol for for different different times times (see text for details). Exposure time to organic contaminants increases from top (blank—no contact) to bottom. bottom. (see text for details). Exposure time to organic contaminants increases from top (blank—no contact) to Dashed Dashed lines indicate the position of diffraction lines typical for bi-hydrated smectites (14.9 and 3.05 Å). lines indicate the position of diffraction lines typical for bi-hydrated smectites (14.9 and 3.05 Å). In any case, the XRD patterns for all samples collected during the sorption isotherms 3.3. Simulations are essentially similar, showing no significant evolution with increasing contact time with The simulation carbamazepine box is shown or paracetamol. in Figure This 3, whereas remarkable Figure stability 4 shows of the the density layer-to-layer profiles distance in- (DPs) calculated along the box z-coordinates for both carbamazepine and dicates that neither paracetamol nor carbamazepine molecules enter the Sap interlayer paracetamol. For both space, PCs, we with consistent observe the two low planes sorptionof capacity water molecules at 2.60–2.65 determined Å andexperiments. from kinetic 6.20–6.30 Å above/below saponite surfaces; a third plane, less pronounced, appears Rather, Ca2+ and Na+ cations persist in these interlayers to compensate for the layerat ~9.15–9.30 Å (Table 1). The water density is constant in the center of the box and corresponds charge deficit, thus keeping constant the overall hydration state of Sap at a given RH. to that of 3.3. Simulations The simulation box is shown in Figure 3, whereas Figure 4 shows the density profiles (DPs) calculated along the box z-coordinates for both carbamazepine and paracetamol. Minerals 2021, 11, 62 For both PCs, we observe two planes of water molecules at 2.60–2.65 Å and 6.20–6.30 7 ofÅ12 above/below saponite surfaces; a third plane, less pronounced, appears at ~9.15–9.30 Å (Table 1). The water density is constant in the center of the box and corresponds to that of liquid liquidwater. water.Sodium Sodiumcations cationsarearealso alsoorganized organizedas astwo twoplanes, planes,with withaawell-defined well-definedfirstfirst atomic atomic plane presenting a higher peak near the upper layer compared to the the plane presenting a higher peak near the upper layer compared to lower lower one, one, with with ~37.7% ~37.7% and ~24.5% and ~24.5% of theoftotal the total number number of Naof Na+ cations, + cations, respectively. respectively. This difference This difference arises arises from a larger number of substitutions in the upper layer (11 out of 19 from a larger number of substitutions in the upper layer (11 out of 19 substitutions), substitutions), thus thus attracting attracting more more cations. cations. TheThe second second cationic cationic planes planes are are at 4.15–4.25 at 4.15–4.25 Å and Å and eacheach con- contains tains ~7% ~7% oftotal of the the total number number of cations. of cations. Overall, Overall, ~75%~75% of charge-compensating of charge-compensating cations cations are are present in these two sets of present in these two sets of planes. planes. Minerals 2021, 11, x FOR PEER REVIEW 8 of 13 Al. Additionally, the carboxyl oxygen coordinates a cation and one aromatic ring with the surface, maximizing van der Waals interactions with the surface. In this mode, N and O atoms are at ~2.0 and 3.2 Å from the layer surface, respectively. The second adsorption mode (λ-mode; Figure 5c) involves only van der Waals interactions between the aromatic rings and the layer surface, with both N and O atoms being at ~7 Å from the layer surface. The two adsorption modes are observed on both layer surfaces, and carbamazepine can thus readily rotate to change its adsorption conformation. In fact, the γ adsorption mode is more frequently observed than the λ-mode on the upper surface (36% and 26% for γ- and λ-modes, respectively) while both modes are observed at the same frequency on the lower surface. This contrast is most likely due to the larger number of substitutions on the Figure Figure3.surface upper 3.Simulation thatbox Simulation boxwith withcarbamazepine promotes at at the the center center as the initial structure cation coordination. carbamazepine structure for for the the production. produc- tion. The paracetamol molecule remains adsorbed onto the upper clay surface during the whole production time. Both the oxygen of the hydroxyl group and the nitrogen of the amine group, in para position, are located at 3.2 Å from the clay surface. This indicates that the paracetamol remains coplanar to the saponite surface, interacting mainly through van der Waals forces and the aromatic ring (see Figure 5a). Interactions with cations through the carboxyl or hydroxyl oxygen atoms of paracetamol are negligible, as indi- cated by the small number of interacting cations (Figure S1 in the Supplementary Materi- als). Carbamazepine first adsorbs to the lower layer surface, then desorbs, and finally adsorbs to the upper surface. Carbamazepine thus spent ~34% of the production time bound to the lower layer surface before desorption, transfer, and adsorption to the upper layer sur- face, where it remained until the end of the production. On both surfaces, there are two peaks for both the amine nitrogen and the carboxyl oxygen atoms. For nitrogen, the first peak is located at 1.95–2.00 Å and the second peak is at 6.65–7.20 Å from the surface oxy- gen; for oxygen, the two peaks are located at 3.15–3.20 Å and 6.85–7.05 Å. For each atom, these two peaks correspond to two adsorption modes (Figure 5b,c). In the γ-mode (Figure Figure Figure 4. 4. Density Density profiles profiles (DPs) (DPs) along along the the z-coordinate z-coordinate (interlayer (interlayer separation orthogonal to basal surfaces). Surface oxygen 5b), carbamazepine interacts with separation the surfaceorthogonal throughtothebasal surfaces). amine groupSurface andoxygen hydrogen (Os) DPs are in red, water oxygen (Os) DPs are in red, water oxygen (Ow) (Ow) in light in light blue, and sodium blue, and sodium cations cations in purple; in purple; for for carbamazepine carbamazepine (left), the (left),tothe snitrogen nitrogens bonds involving negatively charged surface oxygen atoms directly bind the are in blue and the oxygen is in green; for paracetamol (right), the hydroxyl oxygen is in blue and the nitrogen is in green. substituted are in blue and the oxygen is in green; for paracetamol (right), the hydroxyl oxygen is in blue and the nitrogen is in green. Table 1. Density profile data (see Figure 4). Peak coordinates (along z) and densities are given in Å. In parentheses, the percentage of the total atom/molecule density is reported. Paracetamol 1st peak 2nd peak 3rd peak Lower Surf. Upper Surf. Lower Surf. Upper Surf. Lower Surf. Upper Surf. Ow 2.65 (7.5%) 2.65 (7.6%) 6.20 (7.1%) 6.25 (6.9%) 9.15 (6.6%) 9.25 (5.0%) Na+ 2.10 (24.5%) 2.10 (37.1%) 4.25 (7.6%) 4.20 (6.5%) NH 3.2 (97.8%) Minerals 2021, 11, 62 8 of 12 Table 1. Density profile data (see Figure 4). Peak coordinates (along z) and densities are given in Å. In parentheses, the percentage of the total atom/molecule density is reported. Paracetamol 1st peak 2nd peak 3rd peak Lower Surf. Upper Surf. Lower Surf. Upper Surf. Lower Surf. Upper Surf. Ow 2.65 (7.5%) 2.65 (7.6%) 6.20 (7.1%) 6.25 (6.9%) 9.15 (6.6%) 9.25 (5.0%) Na+ 2.10 (24.5%) 2.10 (37.1%) 4.25 (7.6%) 4.20 (6.5%) NH 3.2 (97.8%) OH 3.2 (97.1%) Carbamazepine 1st peak 2nd peak 3rd peak Lower Surf. Upper Surf. Lower Surf. Upper Surf. Lower Surf. Upper Surf. Ow 2.60 (7.3%) 2.65 (7.7%) 6.20 (6.8%) 6.30 (7.0%) 9.15 (8.0%) 9.30 (7.7%) Na+ 2.05 (24.4%) 2.10 (38.3%) 4.20 (6.8%) 4.15 (6.8%) NH2 1.95 (15.7%) 2.00 (36.4%) 7.20 (18.7%) 6.65 (26.0%) O 3.15 (16.2%) 3.20 (36.7%) 7.05 (18.4%) 6.85 (25.7%) The paracetamol molecule remains adsorbed onto the upper clay surface during the whole production time. Both the oxygen of the hydroxyl group and the nitrogen of the amine group, in para position, are located at 3.2 Å from the clay surface. This indicates that the paracetamol remains coplanar to the saponite surface, interacting mainly through van der Waals forces and the aromatic ring (see Figure 5a). Interactions with cations through the carboxyl or hydroxyl oxygen atoms of paracetamol are negligible, as indicated by the small number of interacting cations (Figure S1 in the Supplementary Materials). nerals 2021, 11, x FOR PEER REVIEW 9 of 13 Figure 5. Representative snapshots of the adsorption of paracetamol (a) and carbamazepine adopting the γ-mode (b) and the λ-mode (c) onFigure 5. Representative the clay snapshots surfaces. Distances are in of Å.the adsorption Oxygen atoms of areparacetamol (a)carbons shown in red, and carbamazepine adopt-in light in grey, hydrogens ing the γ-mode (b) and the λ-mode (c) on the clay surfaces. Distances are in Å. Oxygen atoms are grey, silicium in yellow, and sodium cations in pink. Water molecules have been omitted for the sake of clarity. shown in red, carbons in grey, hydrogens in light grey, silicium in yellow, and sodium cations in pink. Water molecules have been omitted for the sake of clarity. Carbamazepine first adsorbs to the lower layer surface, then desorbs, and finally adsorbs to the upper surface. Carbamazepine thus spent ~34% of the production time Figure 6 shows the position of the molecules’ center when adsorbed on the upper bound to the lower layer surface before desorption, transfer, and adsorption to the upper layer surface. Paracetamol remains in diffuse regions free of substitutions, consistent with layer surface, where it remained until the end of the production. On both surfaces, there the lack of interactions with cations. On the contrary, high-density positions of carbamaz- are two peaks for both the amine nitrogen and the carboxyl oxygen atoms. For nitrogen, epine at thethesurface are more first peak localized, is located although at 1.95–2.00 Å carbamazepine and the seconddensity peak isappears lowerÅ from the at 6.65–7.20 than that ofsurface paracetamol (carbamazepine adsorption on the upper surface represents oxygen; for oxygen, the two peaks are located at 3.15–3.20 Å and 6.85–7.05 only Å. For 2/3 of the production each atom, these two peaks correspond to two adsorption modes (Figure to time). The high-density regions are systematically located close 5b,c). In the substitutedγ-mode tetrahedral sites, (Figure where 5b), carbamazepine carbamazepine is retained interacts through with the surfacecation coordina- through the amine group tion, corresponding and hydrogento the bonds adsorption in thenegatively involving γ-mode. When adsorbed charged surfaceinoxygen the λ-mode, atomscar- directly bind bamazepine is more mobile, interacting only through van der Waals interactions with to the substituted Al. Additionally, the carboxyl oxygen coordinates a cation and one the layer surface. In contrast, aromatic for the ring with theγ-mode, surface,carbamazepine maximizing van remains in theinteractions der Waals same adsorption with the surface. site until a cation exchange occurs, which drags the molecule, leading In this mode, N and O atoms are at ~2.0 and 3.2 Å from the layer to the disruption surface, of respectively. the H-bondThe with surface second oxygens.mode (λ-mode; Figure 5c) involves only van der Waals interactions adsorption Minerals 2021, 11, 62 9 of 12 between the aromatic rings and the layer surface, with both N and O atoms being at ~7 Å from the layer surface. The two adsorption modes are observed on both layer surfaces, and carbamazepine can thus readily rotate to change its adsorption conformation. In fact, the γ Figure 5. Representative snapshots of the adsorption of paracetamol (a) and carbamazepine adopt- adsorption mode is more frequently observed than the λ-mode on the upper surface (36% ing the γ-mode (b) and the λ-mode (c) on the clay surfaces. Distances are in Å. Oxygen atoms are and 26% shown in red, carbons for hydrogens in grey, γ- and λ-modes, respectively) in light grey, silicium inwhile yellow,both modes cations and sodium are observed in at the same frequency pink. Water molecules on the have been lowerforsurface. omitted the sakeThis contrast is most likely due to the larger number of of clarity. substitutions on the upper surface that promotes cation coordination. Figure 6 showsFigure 6 shows the position of the the position molecules’of the centermolecules’ center when when adsorbed on theadsorbed upper on the upper layer surface. Paracetamol remains in diffuse regions layer surface. Paracetamol remains in diffuse regions free of substitutions, consistent withfree of substitutions, consistent with the lack of interactions with cations. On the contrary, the lack of interactions with cations. On the contrary, high-density positions of carbamaz- high-density positions of carbamazepine at the surface are more localized, although epine at the surface are more localized, although carbamazepine density appears lower carbamazepine density appears lower than(carbamazepine than that of paracetamol that of paracetamol (carbamazepine adsorption on the upper adsorption surface on the upper represents surface represents only only 2/3 2/3 of the production of the time). Theproduction high-density time). Theare regions high-density systematicallyregions are close located systematically to located close to substituted tetrahedral sites, where carbamazepine substituted tetrahedral sites, where carbamazepine is retained through cation coordina- is retained through cation coordination, corresponding to the adsorption in tion, corresponding to the adsorption in the γ-mode. When adsorbed in the λ-mode, car-the γ-mode. When adsorbed in the λ- mode, carbamazepine is more mobile, interacting only bamazepine is more mobile, interacting only through van der Waals interactions with the through van der Waals interactions layer surface. Inwith the layer contrast, surface. for the γ-mode, In contrast, for theremains carbamazepine γ-mode,incarbamazepine remains in the same the same adsorption site until a cation exchange occurs, which drags the molecule, leading to the disruption of leading to the adsorption site until a cation exchange occurs, which drags the molecule, disruption of the H-bond with surface oxygens. the H-bond with surface oxygens. Figure 6. Saponite upper basal surface and density plots (in blue) of paracetamol (a) and carbamazepine (b) representing Figure 6. Saponite upper basal surface and density plots (in blue) of paracetamol (a) and carbam- their position when adsorbed at the surface. Oxygen atoms are shown in red, silicium in yellow, and aluminum (tetrahedral azepine (b) representing their position when adsorbed at the surface. Oxygen atoms are shown in substitutions) red,insilicium purple.in yellow, and aluminum (tetrahedral substitutions) in purple. In addition, the desorption energies of the two PCs were computed from a free energy In addition, the desorption energies of the two PCs were computed from a free en- profile constructed using an umbrella sampling method to assess the strength of their ergy profile constructed using an umbrella sampling method to assess the strength of their adsorption to the clay surface. PMF curves (Figure S2 in the Supplementary Materials) adsorption to the clay surface. PMF curves (Figure S2 in the Supplementary Materials) show minimal distances between the centers of the molecules and the layer surfaces of 2.95 show minimal distances between the centers of the molecules and the layer surfaces of and 3.12 Å for paracetamol and carbamazepine, respectively. In addition, the desorption 2.95 and 3.12 Å for paracetamol and carbamazepine, respectively. In addition, the desorp- energies were evaluated at 9.6 and 45.2 kJ mol−1 for paracetamol and carbamazepine, tion energies were evaluated at 9.6 and 45.2 kJ mol−1 for paracetamol and carbamazepine, respectively. Although these values must be considered cautiously, they provide a strong respectively. Although these values must be considered cautiously, they provide a strong indication that carbamazepine adsorption may be stronger than that of paracetamol, owing to the presence of cation coordination and H-bond interactions in the former compared to the latter. However, this is not observable from the experimental data obtained in the present study. 4. Conclusions The adsorption of the pharmaceutical compounds carbamazepine and paracetamol onto smectite clays has been investigated using synthetic saponite. The study combines Minerals 2021, 11, 62 10 of 12 kinetic experiments, XRD measurements, and MD simulations with the aim to obtain a complete picture of adsorption events. According to kinetic studies, PC adsorption is somehow limited for both compounds. XRD data show that neither carbamazepine nor paracetamol enter the interlayer space of saponite, suggesting that adsorption occurs only on clay external basal surfaces, leading to limited adsorption. MD simulations investigated the adsorption mechanism on clay basal surfaces, showing interaction resulting mainly from van der Waals interactions between the aromatic rings and the clay surface. For paracetamol, molecule adsorbed parallel to the saponite surface was the only configuration identified. For carbamazepine, however, two adsorption modes were identified. The first one relies also on van der Waals interactions between the aromatic rings and the clay surface (λ-mode), whereas the other adsorption mode (γ-mode) involves stronger interactions. These interactions occur through both coordination to interlayer cations, anchored at substituted tetrahedral sites, and hydrogen bonds with the more acidic oxygen atoms coordinating tetrahedral Al. Supplementary Materials: The following are available online at https://www.mdpi.com/2075-163 X/11/1/62/s1, Figure S1: Radial distribution functions for (a) hydroxyl oxygen atom of paracetamol and sodium cations, (b) carboxyl oxygen atom of paracetamol and sodium cations, and (c) carboxyl oxygen of carbamazepine and sodium cations. In red, the distance at which the rdf’s maximum is found, and in blue, the coordination number computed over the whole simulation; Figure S2: Potential of mean force (kJ mol−1 ) computed for paracetamol (purple curve) and carbamazepine (green curve). Author Contributions: Conceptualization, A.R. and P.M.; methodology, G.C., E.V., B.L. and P.M.; software, G.C.; validation, G.C., E.V., B.L., A.R. and P.M.; formal analysis, E.V., B.L. and P.M.; investigation, G.C., E.V., B.L., A.R. and P.M.; resources, E.V., B.L. and P.M.; data curation, G.C., E.V., B.L., A.R. and P.M.; writing—original draft preparation, P.M.; writing—review and editing, E.V., B.L., and A.R.; visualization, G.C., E.V., B.L., A.R. and P.M.; supervision, P.M.; project administration, A.R. and P.M.; funding acquisition, E.V., B.L., A.R. and P.M. All authors have read and agreed to the published version of the manuscript. Funding: Financial support from the “PEPS Lyon CNRS 2016: Risques et environnement, project ModExSol” is acknowledged. MINECO (project CTQ2017-89132-P) and DIUE (project 2017SGR1323) are acknowledged for financial support. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data is contained within the article or supplementary material. Acknowledgments: Calculations were performed using HPC resources from GENCI- [IDRIS: gen7662]. A.R. is indebted to the “Ramón y Cajal” program. Conflicts of Interest: The authors declare no conflict of interest. References 1. Clarke, B.O.; Smith, S.R. Review of ‘emerging’ organic contaminants in biosolids and assessment of international research priorities for the agricultural use of biosolids. Environ. Int. 2011, 37, 226–247. [CrossRef] [PubMed] 2. Lapworth, D.J.; Baran, N.; Stuart, M.E.; Ward, R.S. Emerging organic contaminants in groundwater: A review of sources, fate and occurrence. Environ. Pollut. 2012, 163, 287–303. [CrossRef] [PubMed] 3. Vodyanitskii, Y.N.; Yakovlev, A.S. Contamination of soils and groundwater with new organic micropollutants: A review. Eurasian Soil Sci. 2016, 49, 560–569. [CrossRef] 4. Heberer, T. Occurrence, fate, and removal of pharmaceutical residues in the aquatic environment: A review of recent research data. Toxicol. Lett. 2002, 131, 5–17. [CrossRef] 5. Ternes, T.A. Occurrence of drugs in German sewage treatment plants and rivers1Dedicated to Professor Dr. Klaus Haberer on the occasion of his 70th birthday.1. Water Res. 1998, 32, 3245–3260. [CrossRef] 6. Wilken, R.D.; Ternes, T.A.; Heberer, T. Pharmaceuticals in Sewage, Surface and Drinking Water in Germany. In Security of Public Water Supplies; Deininger, R.A., Literathy, P., Bartram, J., Eds.; Springer: Dordrecht, The Netherlands, 2000; pp. 227–240. ISBN 978-0-7923-6122-0. Minerals 2021, 11, 62 11 of 12 7. Castiglioni, S.; Bagnati, R.; Fanelli, R.; Pomati, F.; Calamari, D.; Zuccato, E. Removal of Pharmaceuticals in Sewage Treatment Plants in Italy. Environ. Sci. Technol. 2006, 40, 357–363. [CrossRef] 8. Loraine, G.A.; Pettigrove, M.E. Seasonal Variations in Concentrations of Pharmaceuticals and Personal Care Products in Drinking Water and Reclaimed Wastewater in Southern California. Environ. Sci. Technol. 2006, 40, 687–695. [CrossRef] 9. Reemtsma, T.; Weiss, S.; Mueller, J.; Petrovic, M.; González, S.; Barcelo, D.; Ventura, F.; Knepper, T.P. Polar Pollutants Entry into the Water Cycle by Municipal Wastewater: A European Perspective. Environ. Sci. Technol. 2006, 40, 5451–5458. [CrossRef] 10. Al-Rifai, J.H.; Gabelish, C.L.; Schäfer, A.I. Occurrence of pharmaceutically active and non-steroidal estrogenic compounds in three different wastewater recycling schemes in Australia. Chemosphere 2007, 69, 803–815. [CrossRef] 11. Ellis, J.B. Pharmaceutical and personal care products (PPCPs) in urban receiving waters. Environ. Pollut. 2006, 144, 184–189. [CrossRef] 12. Yu, J.T.; Bouwer, E.J.; Coelhan, M. Occurrence and biodegradability studies of selected pharmaceuticals and personal care products in sewage effluent. Agric. Water Manag. 2006, 86, 72–80. [CrossRef] 13. Grossberger, A.; Hadar, Y.; Borch, T.; Chefetz, B. Biodegradability of pharmaceutical compounds in agricultural soils irrigated with treated wastewater. Environ. Pollut. 2014, 185, 168–177. [CrossRef] [PubMed] 14. Lin, K.; Gan, J. Sorption and degradation of wastewater-associated non-steroidal anti-inflammatory drugs and antibiotics in soils. Chemosphere 2011, 83, 240–246. [CrossRef] [PubMed] 15. Kolpin, D.W.; Furlong, E.T.; Meyer, M.T.; Thurman, E.M.; Zaugg, S.D.; Barber, L.B.; Buxton, H.T. Pharmaceuticals, Hormones, and Other Organic Wastewater Contaminants in U.S. Streams, 1999−2000: A National Reconnaissance. Environ. Sci. Technol. 2002, 36, 1202–1211. [CrossRef] [PubMed] 16. Snyder, S.A.; Adham, S.; Redding, A.M.; Cannon, F.S.; DeCarolis, J.; Oppenheimer, J.; Wert, E.C.; Yoon, Y. Role of membranes and activated carbon in the removal of endocrine disruptors and pharmaceuticals. Desalination 2007, 202, 156–181. [CrossRef] 17. Chefetz, B.; Mualem, T.; Ben-Ari, J. Sorption and mobility of pharmaceutical compounds in soil irrigated with reclaimed wastewater. Chemosphere 2008, 73, 1335–1343. [CrossRef] 18. Salvia, M.V.; Experton, J.; Geandel, C.; Cren-Olivé, C.; Vulliet, E. Fate of pharmaceutical compounds and steroid hormones in soil: Study of transfer and degradation in soil columns. Environ. Sci. Pollut. Res. 2014, 21, 10525–10535. [CrossRef] 19. Salvia, M.V.; Cren-Olivé, C.; Vulliet, E. Statistical evaluation of the influence of soil properties on recoveries and matrix effects during the analysis of pharmaceutical compounds and steroids by quick, easy, cheap, effective, rugged and safe extraction followed by liquid chromatography–tandem mass spectrometry. J. Chromatogr. A 2013, 1315, 53–60. [CrossRef] 20. Stathi, P.; Litina, K.; Gournis, D.; Giannopoulos, T.S.; Deligiannakis, Y. Physicochemical study of novel organoclays as heavy metal ion adsorbents for environmental remediation. J. Colloid Interface Sci. 2007, 316, 298–309. [CrossRef] 21. Abollino, O.; Aceto, M.; Malandrino, M.; Sarzanini, C.; Mentasti, E. Adsorption of heavy metals on Na-montmorillonite. Effect of pH and organic substances. Water Res. 2003, 37, 1619–1627. [CrossRef] 22. Chang, P.H.; Li, Z.; Jiang, W.T.; Sarkar, B. Clay minerals for pharmaceutical wastewater treatment. In Modified Clay and Zeolite Nanocomposite Materials; Elsevier: Amsterdam, The Netherlands, 2019; pp. 167–196. ISBN 978-0-12-814617-0. 23. Thiebault, T.; Boussafir, M. Adsorption Mechanisms of Psychoactive Drugs onto Montmorillonite. Colloid Interface Sci. Commun. 2019, 30, 100183. [CrossRef] 24. Thiebault, T.; Boussafir, M.; Le Forestier, L.; Le Milbeau, C.; Monnin, L.; Guégan, R. Competitive adsorption of a pool of pharmaceuticals onto a raw clay mineral. RSC Adv. 2016, 6, 65257–65265. [CrossRef] 25. Mignon, P.; Ugliengo, P.; Sodupe, M. Theoretical Study of the Adsorption of RNA/DNA Bases on the External Surfaces of Na+ -Montmorillonite. J. Phys. Chem. C 2009, 113, 13741–13749. [CrossRef] 26. Mignon, P.; Sodupe, M. Structural Behaviors of Cytosine into the Hydrated Interlayer of Na+ -Montmorillonite Clay. An ab Initio Molecular Dynamics Study. J. Phys. Chem. C 2013, 117, 26179–26189. [CrossRef] 27. Mignon, P.; Navarro-Ruiz, J.; Rimola, A.; Sodupe, M. Nucleobase Stacking at Clay Edges, a Favorable Interaction for RNA/DNA Oligomerization. ACS Earth Space Chem. 2019, 3, 1023–1033. [CrossRef] 28. Mignon, P.; Corbin, G.; Le Crom, S.; Marry, V.; Hao, J.; Daniel, I. Adsorption of nucleotides on clay surfaces: Effects of mineral composition, pH and solution salts. Appl. Clay Sci. 2020, 190, 105544. [CrossRef] 29. Galicia-Andrés, E.; Tunega, D.; Gerzabek, M.H.; Oostenbrink, C. On glyphosate-kaolinite surface interactions. A molecular dynamic study. Eur. J. Soil Sci. 2020. [CrossRef] 30. Aristilde, L.; Marichal, C.; Miéhé-Brendlé, J.; Lanson, B.; Charlet, L. Interactions of Oxytetracycline with a Smectite Clay: A Spectroscopic Study with Molecular Simulations. Environ. Sci. Technol. 2010, 44, 7839–7845. [CrossRef] 31. Hamilton, D.L.; Henderson, C.M.B. The preparation of silicate compositions by a gelling method. Mineral. Mag. J. Mineral. Soc. 1968, 36, 832–838. [CrossRef] 32. Robert, J.L.; Beny, J.M.; Ventura, G.D.; Hardy, M. Fluorine in micas: Crystal-chemical control of the OH-F distribution between trioctahedral and dioctahedral sites. Eur. J. Mineral. 1993, 5, 7–18. [CrossRef] 33. Bergaoui, L.; Lambert, J.F.; Franck, R.; Suquet, H.; Robert, J.L. Al-pillared saponites. Part 3.—Effect of parent clay layer charge on the intercalation–pillaring mechanism and structural properties. J. Chem. Soc. Faraday Trans. 1995, 91, 2229–2239. [CrossRef] 34. Michot, L.; Villieras, F. Assessment of surface energetic heterogeneity of synthetic Na-saponites. The role of layer charge. Clay Miner. 2002, 37, 39–57. [CrossRef] Minerals 2021, 11, 62 12 of 12 35. Loos, R.; Gawlik, B.M.; Locoro, G.; Rimaviciute, E.; Contini, S.; Bidoglio, G. EU-wide survey of polar organic persistent pollutants in European river waters. Environ. Pollut. 2009, 157, 561–568. [CrossRef] [PubMed] 36. Loos, R.; Locoro, G.; Comero, S.; Contini, S.; Schwesig, D.; Werres, F.; Balsaa, P.; Gans, O.; Weiss, S.; Blaha, L.; et al. Pan-European survey on the occurrence of selected polar organic persistent pollutants in ground water. Water Res. 2010, 44, 4115–4126. [CrossRef] 37. López-Pacheco, I.Y.; Silva-Núñez, A.; Salinas-Salazar, C.; Arévalo-Gallegos, A.; Lizarazo-Holguin, L.A.; Barceló, D.; Iqbal, H.M.N.; Parra-Saldívar, R. Anthropogenic contaminants of high concern: Existence in water resources and their adverse effects. Sci. Total Environ. 2019, 690, 1068–1088. [CrossRef] 38. Woods, R.J.; Chappelle, R. Restrained electrostatic potential atomic partial charges for condensed-phase simulations of carbohy- drates. J. Mol. Struct. THEOCHEM 2000, 527, 149–156. [CrossRef] 39. Cygan, R.T.; Liang, J.J.; Kalinichev, A.G. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. J. Phys. Chem. B 2004, 108, 1255–1266. [CrossRef] 40. Berendsen, H.J.C.; Postma, J.P.M.; Van Gunsteren, W.F.; Hermans, J. Interaction models for water in relation to protein hydration. Jerusalem Symp. Quantum Chem. Biochem. 1981, 14, 331–342. 41. Dang, L.X. Mechanism and Thermodynamics of Ion Selectivity in Aqueous Solutions of 18-Crown-6 Ether: A Molecular Dynamics Study. J. Am. Chem. Soc. 1995, 117, 6954–6960. [CrossRef] 42. Berendsen, H.J.C.; Grigera, J.R.; Straatsma, T.P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269–6271. [CrossRef] 43. Wang, J.; Wolf, R.; Caldwell, J.; Kollman, P.; Case, D. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174. [CrossRef] [PubMed] 44. Thyveetil, M.A.; Coveney, P.V.; Greenwell, H.C.; Suter, J.L. Computer simulation study of the structural stability and materials properties of DNA-intercalated layered double hydroxides. J. Am. Chem. Soc. 2008, 130, 4742–4756. [CrossRef] [PubMed] 45. Swadling, J.B.; Suter, J.L.; Greenwell, H.C.; Coveney, P.V. Influence of Surface Chemistry and Charge on Mineral–RNA Interactions. Langmuir 2013, 29, 1573–1583. [CrossRef] [PubMed] 46. Swadling, J.B.; Coveney, P.V.; Greenwell, H.C. Clay Minerals Mediate Folding and Regioselective Interactions of RNA: A Large-Scale Atomistic Simulation Study. J. Am. Chem. Soc. 2010, 132, 13750–13764. [CrossRef] [PubMed] 47. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. [CrossRef] 48. Ryckaert, J.P.; Ciccotti, G.; Berendsen, H.J.C. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327–341. [CrossRef] 49. 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] 50. 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. Comput. Phys. Commun. 2009, 180, 1961–1972. [CrossRef] 51. The PLUMED consortium Promoting transparency and reproducibility in enhanced molecular simulations. Nat. Meth. 2019, 16, 670–673. [CrossRef] 52. Ferrage, E.; Lanson, B.; Malikova, N.; Plançon, A.; Sakharov, B.A.; Drits, V.A. New Insights on the Distribution of Interlayer Water in Bi-Hydrated Smectite from X-ray Diffraction Profile Modeling of 00 l Reflections. Chem. Mater. 2005, 17, 3499–3512. [CrossRef] 53. Dazas, B.; Lanson, B.; Delville, A.; Robert, J.L.; Komarneni, S.; Michot, L.J.; Ferrage, E. Influence of Tetrahedral Layer Charge on the Organization of Interlayer Water and Ions in Synthetic Na-Saturated Smectites. J. Phys. Chem. C 2015, 119, 4158–4172. [CrossRef]