Long-Range Ionic and Short-Range Hydration Effects Govern Strongly Anisotropic Clay Nanoparticle Interactions

The aggregation of clay particles in aqueous solution is a ubiquitous everyday process of broad environmental and technological importance. However, it is poorly understood at the all-important atomistic level since it depends on a complex and dynamic interplay of solvent-mediated electrostatic, hydrogen bonding, and dispersion interactions. With this in mind, we have performed an extensive set of classical molecular dynamics simulations (included enhanced sampling simulations) on the interactions between model kaolinite nanoparticles in pure and salty water. Our simulations reveal highly anisotropic behavior, in which the interaction between the nanoparticles varies from attractive to repulsive depending on the relative orientation of the nanoparticles. Detailed analysis reveals that at large separation (>1.5 nm), this interaction is dominated by electrostatic effects, whereas at smaller separations, the nature of the water hydration structure becomes critical. This study highlights an incredible richness in how clay nanoparticles interact, which should be accounted for in, for example, coarse-grained models of clay nanoparticle aggregation.

The bulk structure of kaolinite is known and well characterised. For the present work, we took as reference the crystallographic structure of Bish [1] [2]. The relaxed structures obtained using the chosen ClayFF force field or using DFT-based approaches are in close agreement with the crystallographic experimental bulk structure of kaolinite. The crystallographic structure of kaolinite nanoparticles is instead not immediately available. However, from the literature there are two very stable basal surfaces, and the edge surfaces correspond to the following planes: (0 1 0), (0 -1 0), (1 1 0), (-1 -1 0), (1 -1 0), and (-1 1 0). Cleavage along these planes has also been studied in the literature, using force fields and DFT, and precise indications on the cleaving and hydrogenation of dangling bonds are provided, for instance, in [3] and references therein, resulting in five-fold coordination of edge Al atoms. Here we followed those prescriptions, and considered a 3-layer (∼ 20Å) hexagonal kaolinite nanoparticle with side lengths of ∼ 18Å. The nanoparticle was created in the the VESTA package [4] and dangling bonds on edges were hydrogenated using a Python script written explicitly to fulfil the indications given in the literature. The surface area for each edge is roughly A edge ∼ 18 × 20 = 360Å 2 , and the surface area of the basal surfaces is roughly The ratio between the face and the edge surface of the simulated particles is A edge A edge ∼ 2.3. We expect that face-to-face interfaces are, in a first approximation, proportional to A face and the face-to-edge interfaces are proportional to A edge .
The results shown in the manuscript are obtained from classical molecular dynamics simulations employing semi-empirical force fields, are discussed in the computational details. In each simulation cell there are two kaolinite nanoparticles and 10,000 to 11,000 water molecules. In line with ref. [5], the simulation cell is roughly 60 × 60 × 100Å 3 for PMF calculations of face-to-face orientation, and 60 × 60 × 110Å 3 for PMFs calculations of face-to-edge orientation. Periodic boundary conditions are employed in all three dimensions. An image of the nanoparticles in face to face configuration is shown in Figure S1. S4 Confidential FIG. S1. A representative simulation snapshot of the face-to-face case in the gibbsite-siloxane orientation. The blue box is the simulation cell, where periodic boundary conditions are employed. Inside the box, two nanoparticles are surrounded by 10,000 water molecules, coloured in cyan. Other colors are the same as those shown in Figure 1 in the manuscript. The distance among the two nanoparticles is evaluated during the MD simulation by measuring the projection on the z-axis of the distance among the center of mass of the heavy atoms in the two nanoparticles, which is indicated as ∆z com in the image. The umbrella sampling constraints in our simulations are enforced on ∆z com , and in this specific orientation we performed MD simulations with ∆z com ranging from 2 nm to 5 nm. The PMF curves have been plotted as functions of ∆Z, which is the minimum distance among the heavy atoms in the two particles, because it facilitates the comparison among the curves of different orientations. Indeed, ∆z com and ∆Z differ by an offset, which depends on the specific orientation among the nanoparticles.

S2. DENSITY FUNCTIONAL THEORY TESTS
In order to validate the reliability of the force field parametrization employed in the present work (described in the methods section of the manuscript), we performed a set of simulations using ab-initio approaches. Among the latter methods, the most affordable is density functional theory (DFT). DFT is able to provide evaluations in good agreement with experiment or higher level of theories, as quamtum Monte Carlo, in systems including clay, water and ions, as long as it employs a correction for the Van der Waals (vdW) [6]. There are many different flavours of DFT and of associated vdW-correction schemes, and some setups proves more reliable than others. Based on some benchmarks reported in [6], we selected a good setup for DFT and we performed DFT- Adsorption was examined on a single layer of kaolinite, a system with 2D periodicity along the A and B axes. The simulated supercell was 1 × 2 the conventional unit cell of bulk kaolinite (ca. 10.38Å× 9.01Å) and is comprised of 8 aluminium, 8 silicons, 36 oxygen, and 16 hydrogen S6 atoms for the kaolinite slab, plus the atoms of the adsorbed molecule or ion. We used a C axis of 22.25Å, yielding ca. 15Å of vacuum in the direction orthogonal to the slab. We relaxed the atomic configurations, keeping the cell axes fixed, in all the systems under consideration: the free standing slab, and the slab plus the adsorbate on either the gibbsite or the siloxane face, being the adsorbate either a water molecule, a Na + ion or a Cl − ion. DFT simulations were performed with the CP2K [7] software package using the Quickstep algorithm, the PBE-D3(BJ) [8][9][10] functional, Goedecker-Teter-Hutter pseudopotentials [11,12], a 1000 Ry plane-wave cut-off, a relative cut-off of 60 Ry, TZVP molecularly optimised basis set [13] on H, O, and Si, and DZVP molecularly optimised basis set on Al and Cl (because the TZVP was not available). The plane-wave cutoffs were carefully tested on the adsorption of water on both faces of kaolinite, where the basis set error with the selected basis set is < 1 meV. It was not possible to make a comparable convergence test on the localised basis sets, as some atomic species only the SZV and the DZVP are provided. Thus, we compared the adsorption energy of water on both surfaces with the results from the VASP [14] package, using the same force field, the PAW potentials [15] and converged plane wave cutoffs.
The VASP and CP2K results agree within 1 kcal/mol. In the simulations where the adsorbate is an ion the overall system is charged, either positively or negatively, so a uniform compensating background charge is added, as it is standard practice in DFT.
Each adsorption energy E ads was estimated as the difference between the energy of the optimised configuration E opt and the energy of the systems far away E far : E ads = E opt − E far . In particular, E far is the energy of the system of the relaxed structure of the free standing slab plus the adsorbate in the middle of the vacuum region (i.e., over 7Å away from both the faces) [16]. The E ads computed in this way are reported in Table S1. They are compared with the values obtained employing the classical force field setup that we employed in the rest of the project, namely ClayFF and the SPC/E water model. In the FF approach we used the same simulation cell as in DFT, but we relaxed the atomic positions. The DFT and the FF relaxed structures are in good agreement: they is a qualitatively good matching, and the main differences are little interatomic deviations.
The values of E ads computed using ClayFF+SPC/E seem to always underestimate the interaction strength if compared with DFT, however it yields to right ranking. It is worth to consider how much the DFT results depend on the particular setup that we have used. In this respect, the values from VASP, using the same PBE-D3 functional, gives -17.7 kcal/mol and -5.7 kcal/mol on W@G and W@S, respectively. This shows that the cutoffs are at convergence and there is little dependence on the employed pseudopotentials. The uncertainty due to the choice of the functional and the vdW correction is much larger, as shown for instance in Zen et al. [6]. Ref. [6] also shows TABLE S1. Adsorption energy E ads , in kcal/mol, of a water (W) molecule, a chlorine (Cl) ion or a sodium (Na) ion on either the gibbsite (G) or the siloxane (S) faces. Values have been obtained using the LAMMPS package [17] on the classical force field ClayFF and the SPC/E water model, and using the CP2K package [7] on DFT simulations with the PBE-D3 [8][9][10]   . Protocol that we developed and used to evaluate the PMF of two kaolinite nanoparticles. The first step is to decide the relative orientation of the nanoparticles, and generate the corresponding system. The second step is the equilibration of the system, in the NVT ensemble, for at least 1 ns. The third step is to introduce an external driving force on the upper nanoparticle driving it from far away to very close to the lower nanoparticle. We set this external force in order to move the nanoparticle 1 nm each ns. We are interested in a region of ∆Z ranging from 0 to 3 nm, thus step 3 requires 3 ns of MD simulation. In step 4, a subset of the configurations sampled in the pulling dynamics are chosen and used as starting configurations of n umbrella sampling (US) simulations. Each US simulation samples the phase space on a specific predefined interparticle distance. This is achieved by imposing to ∆z com a harmonic force potential with force constant k i and center ∆ i . In step 5 the n US MD simulations are performed in an NVT ensemble. Finally, the outputs of the n US simulations are used to check the sampling of the ∆Z and build the PMF, using the Weighted Histogram Analysis Method (WHAM). Black boxes in the image correspond to running MD simulations. Further details are given in section S3.
The first step of the scheme was performed using the Packmol package [18] to create the structure of the two nanoparticles in a specific orientation and distance, and of the water molecules around them. Afterwards, we used the Moltemplate package [19] to generate the LAMMPS inputs. The We have checked that they can be performed efficiently using approx. 100 CPU cores. The pulling MD simulation in step 3 is aimed at providing good initial configurations for the following US simulations. The pulling has to be slow enough that the water particles have time to equilibrate around the nanoparticles. A drift velocity of 1 nm/ns was observed to be good. In step 4 a decision on how many US simulations to perform has to be made, which is related to the choice of force constant values k i and the separation (i.e., the value of |∆ i − ∆ i+1 |) among them. We tested that an effective choice is to use a different setup in the regions corresponding to ∆Z larger and smaller of ca. 1 nm. Indeed, the PMF is quite smooth when the nanoparticles are separated by more than 1 nm, and we used a k i of 10 kcal/mol and a separation |∆ i − ∆ i+1 | of around 0.7Å. When the nanoparticles are close the PMF can have steep oscillations, so we used a k i of 50 kcal/mol and a separation |∆ i − ∆ i+1 | of 0.2Å. S10 S4.

TEMPERATURE DEPENDENCE OF PMF
We consider here the dependence on the temperature of the PMF profiles. All the PMFs considered in the previous figures have been obtained at a temperature T=350 K. We performed an additional PMF evaluation at an higher temperature, namely the gibbsite-siloxane orientation at T=450 K. The comparison between the PMFs at T=350 K and T=450 K is shown in figure S5. The comparison indicates that there is not a big sensitivity of the PMF on the temperature, at least for this specific orientation and specific modelling approach employed. The local minima and the barriers are only slightly different at the two temperatures considered.

S5. PMF PER-SURFACE-AREA PROFILES
Since kaolinite nanoparticles can have different sizes from those simulated here, in Figure S6 we have plotted the PMF per surface area. We cannot assume a precise linear dependence between the PMF and the surface area at the interface, because the borders of the interface are not negligible compared to the interface itself. However, in the first instance the assumption of linearity could be considered acceptable. Assuming it, the comparison of the curves in Figure S6 indicates that the "flatness" of the face-to-edge interfaces does not follow only the interface surface area and is an intrinsic property of the nanoparticle interactions.

S7. PARTICLE-PARTICLE INTERACTIONS
In Figure S8 we plot the decomposition of particle-particale interaction into vdW and Coulomb contributions using fix group/group command implemented in LAMMPS [17]. Due to the similar charge distributions of the two approaching surfaces, the results show a strong repulsive Coulomb interaction between the two particles, which is opposite to what has been observed in the case of the gibbsite-siloxane orientation (see Figure 3 of the main text). S8. Results of decomposition of particle-particle interactions: vdW and Coulomb potentials as a function of the distance between the two particles, calculated for siloxane-siloxane (a) and gibbsite-gibbsite (b) orientations, respectively. Noting that the Coulomb potential values are highly positive in both cases due to symmetrical charge distributions on the two approaching particles. S15

S8. SOLVATION FREE ENERGY
To obtain the solvation free energy, we applied the method described in [20][21][22], which is based on the quantification of time-averaged water density fluctuations. These are calculated via the probability of observing N molecules within a 6.6Å diameter spherical volume, v, located at the center of the first hydration film.
Then the solvation free energy is calculated based on the relation: Figure S9 depicts the solvation free energies at various distances. Consistently with the other results presented in this manuscript, we observed a noticeable difference in the solvation free energy between different particle-particle separations in the cases of siloxane-siloxane and gibbsitegibbsites orientations. Whereas, the solvation free energy obtained in the case of gibbsite-siloxane appears to remain constant (see Figure S9c, f, and g). The minimal change in the calculated solvation free energy indicates agglomeration along the gibbsite-siloxane orientation incurs a lower energy cost than along the gibbsite-gibbsite and siloxane-siloxane orientations.
FIG. S9. In-plane solvation free energy distributions obtained for the siloxane and gibbsite surfaces in the siloxane-siloxane (SS), gibbsite-gibbsite (GG), and gibbsite-siloxane (GS) orientated systems. The results were obtained from density fluctuations of water within a 6.6Å spherical volume located in middle of the first hydration layer of the siloxane and gibbsite surfaces (along the X-Y plane). S17 Figure S10 shows the density of water oxygens and hydrogen atoms along the Z direction of the simulation box (perpendicular to the basal planes of the kaolinite particles) at various separation distances, as obtained from the MD simulations in pure water solution. The upper two profiles correspond to the local minima of the PMF. Figure S10(a) for siloxane-siloxane shows that the density of the oxygen atoms is bimodal despite the interface being large enough to accommodate only one layer of water. The density profiles of hydrogen atoms in water molecules in Figure S10(a) exhibits three closely spaced peaks, with the middle peak roughly two times as large as the side peaks.

S9. DENSITY PROFILES OF WATER OXYGEN AND HYDROGEN ATOMS
In contrast, in the case of gibbsite-gibbsite, we observed in Figure S10(b) that the density peaks of water oxygen atoms merge into a single sharp peak and the hydrogen atoms density is again three-modal. Based on the positions of the hydrogen and oxygen peaks, we infer that one of the water molecules' OH groups points toward the two surfaces and one is parallel to the surface.
It is worth noting that the density profiles for siloxane-siloxane and gibbsite-gibbsite are highly symmetrical, as comes from the up-down inversion symmetry of the system.
In the case of gibbsite-siloxane orientation, Figure S10(c) shows asymmetric peaks; the oxygen atom density is unimodal but skewed, and the hydrogen atom density has one little peak close to the gibbsite surface and a big peak close to siloxane. These distributions are qualitatively similar to those of oxygen and hydrogen atoms in the proximity of gibbsite in the large ∆Z case (see Figure   S10h and S10i). This results in a minimal cost required for merging the two hydration films at close inter-particle distances, which is consistent with the results calculated for water orientation (see Figure 4 in the manuscript).
These results indicate that in the cases of siloxane-siloxane and gibbsite-gibbsite orientations, there is a significant change in the relative distribution of water oxygen and hydrogen atoms within the hydration layers as the two particles approach each other. In contrast, we do not observe such a significant change in the case of gibbsite-siloxane orientation. This yields a negligible energy cost for merging the two hydration films at close inter-particle distances for the gibbsite-siloxane orientation, but a penalty for the siloxane-siloxane and gibbsite-gibbsite orientations. This is consistent with the analysis of the water orientation (see Figure 4 in the manuscript and the corresponding discussion).

S10. DENSITY OF WATER
We consider here the average density of water molecules between two approaching particles as a function of particle-particle separation. The results, shown in Figure S11, indicate that when the approaching surfaces contain at least one gibbsite surface (i.e., gibbsite-siloxane and gibbsitegibbsite), the positions of the free energy wells correspond to the distances at which the number of confined water molecules increases. This is almost certainly due to the fact that water molecules prefer to form hydrogen bonds with gibbsite surfaces (see Figure S12). In the siloxane-siloxane orientation, on the other hand, water-enriched regions correspond to energy barriers. This could be because the siloxane surfaces are more hydrophobic, which do not support increased water density within the hydration layer.

Normalized density
FIG. S11. Average density of water molecules between the two approaching surfaces as a function of the inter-particle distance normalised by the density of water at a far distance. The locations of free energy wells and barriers are displayed in green and orange, respectively. Obtained from simulations with pure water. S20 S11. HYDROGEN BONDS In Figure S12

S12. DIFFUSION COEFFICIENT OF WATER
The self-diffusion of water molecules in two different simulation boxes has been estimated and shown as a function of the coordinate Z in Figure S13. In this calculation, the diffusion coefficient values were calculated using Einstein relation in mean square displacement (MSD) curve for water molecules belonging to a specific slab of 20.0 × 20.0 × 3.0Å along the Z direction. Points on the left of the figure correspond to water molecules forming the hydration film of one of the nanoparticles.
At a distance of more than ca. 1.5 nm the diffusion of water is the same as in bulk. In contrast, close to the kaolinite face the diffusion of water is drastically reduced, by a factor 2. , and compared to the results of the bulk water at the same condition. The distances between the two particles in these calculations are 30Å. The results show that diffusion coefficient of water in the middle region between the two particles (around 35-40Å) reaches the bulk value and the effect of the simulation box size is minimal.

S13. PARTICLE DYNAMICS IN SOLUTION
Particle diffusion in solution is one of the factors determining the kinetics of particle aggregation.
Analysis of the umbrella sampling trajectories yields estimates regarding the local diffusion of the kaolinite particles. To extract this information we apply the formalism proposed by Woolf and Roux [23][24][25]: where dz is the average distance along the Z direction between two kaolinite particles under harmonic restraint umbrella sampling simulation at a specific window/distance, Var(dz) = dz 2 − dz 2 is its variance, and τ dz is its correlation time defined as: with δdz(t) = dz(t) − dz .
Using this approach, we extract the diffusion coefficients calculated for the kaolinite particle as a function of the distance between the two particles at different orientations in pure and saline water.