A structural coarse-grained model for clays using simple iterative Boltzmann inversion

The environmental impact of cesium adsorption and diffusion into various types of naturally occurring layered silicates has received renewed interest in recent years, especially in the aftermath of the Fukushima Daiichi Nuclear Power Plant disaster.^1–3 One of the most environmentally threatening products of nuclear fission is cesium-137, both because of its relatively long half-life (30.2 years) and its affinity for mineral surfaces, which prevents it from leaching from surface soils.^4,5 Cesium is strongly and irreversibly adsorbed to various clay surfaces in the presence of other ions and can slowly diffuse into the bulk volume of both anhydrous and hydrated layered silicates.^6–8 Due to its intermediate half-life and its relative abundance as a nuclear decay product, cesium-137 can contaminate environmental sites with dangerous levels of radiation for decades, while most other fallout isotopes may only present a threat on the order of months or years.⁹ 

Cesium diffuses deep into naturally occurring clays and displaces other types of ions normally found in the interlayer, such as potassium, sodium, and calcium. Because of its ability to selectively exchange radiocesium, illite and similar clays have been investigated for the possibility of remediating radioactive plumes of cesium.¹⁰ Despite extensive experimental study of ion adsorption at the frayed clay edge and exchange of ions in the clay interlayer,^7,8,11–13 the exact mechanism of cesium uptake remains elusive.¹⁴ It is especially unclear how cesium displaces potassium within the interlayer far from the edge in anhydrous interlayers of clays such as illite,^11,15–17 despite very large barriers to ion diffusion. Some groups have suggested that only hydrated ions within the interlayer or near a clay edge should be capable of exchange;¹⁸ however, this exchange mechanism is very thermodynamically unfavorable in bulk illite, since both potassium and cesium strongly favor anhydrous interlayers.¹⁹ Moreover, cesium is found to have a strong affinity for the weathered clay edge, further suggesting that interlayer hydration is likely not a necessary step either for cesium binding or penetration into the interlayer.^14,20–24

Recent experimental^7,8 and computational work²⁵ supports the direct exchange by diffusion for large ions in the interlayer of several clay types without repeated hydration and dehydration at the edge sites. Instead, cesium has been hypothesized to first bind to frayed clay edges and then more slowly exchange with the potassium naturally found in the interlayer through a diffusive mechanism.^17,26 Once potassium begins exchanging for cesium at the weathered edge, it has been proposed that there is a thermodynamic and kinetic driving force to displace additional potassium ions, resulting in an accelerated replacement of potassium ions for cesium.^27,28 Despite this proposed feed-forward mechanism, it has been observed that the rates of ion exchange over moderate time scales in neighboring interlayers can vary drastically.^7,8 Although some thermodynamic arguments have been proposed to explain the stability of interstratified clay particles in other clay types,^29–31 the physical reason for the disparity in neighboring interlayer exchange rates in illite is still not entirely clear.

In order to adequately predict the ability of micaceous minerals both to uptake radioactive cesium and its susceptibility to remobilization, it is necessary to understand and model the mechanisms that drive the adsorption and exchange of various ionic species in, for example, anhydrous illite. Recent computational work from our group using the classical ClayFF force field and density functional theory (DFT) has indicated that the presence of cesium ions within potassium interlayers creates mechanical forces that significantly increase the interlayer spacing.²⁵ This expansion in turn results in lower rate coefficients for ion exchange by 6–10 orders of magnitude, which accelerates the diffusion of potassium ions.²⁵ However, the problem of characterizing interstratified particles by definition requires the evaluation of diffusion barriers involving multiple interlayers. Probing the energy landscape of large atomistic simulations of this type is computationally expensive and motivates the development of a coarse-grained (CG) model that is capable of reproducing the energy barriers to ion diffusion near an exchange front and in interstratified clay particles.

Previous CG simulations of clay include both continuous and discrete models of the clay interlayer.^32–35 In one study, Marry et al.³² investigated a montmorillonite system with a hydrated interlayer. Their CG model consisted of two uniformly charged plates representing the clay layers held a fixed distance apart. Water was modeled implicitly, and both cations and anions in solution were modeled continuously as an effective density that was allowed to vary over space. This model accounted for the excluded volume of the ions using the mean spherical approximation, an advance on previous continuous models of hydrated clay systems. Using DFT, the grand potential was minimized to derive the density of ions as a function of distance from the clay surface for sodium ions. These types of CG models can be highly accurate for the calculation of continuous properties such as tensile strength or average interlayer spacing. Several groups have also successfully modeled clay and clay-polymer systems by coarse-graining entire clay layers into one single particle.^36–38 However, continuous, ultra coarse-grained models are fundamentally incapable of capturing site-specific effects such as ion binding and due to mean-field approximations, would not be expected to reproduce the observed variation in the ion binding energy landscape in anhydrous clay.²⁵ 

Other groups have developed more detailed coarse-grained clay models that map multiple atoms to single pseudo-atoms, a common coarse-graining technique.³⁹ In the model of clay-polymer composites developed by the Coveney group, each montmorillonite clay layer is represented by bonded pseudo-atoms corresponding to ion binding sites.^40,41 The pseudo-atoms in the Suter et al. model are bonded with harmonic potentials derived from an iterative Boltzmann inversion (IBI) of the corresponding all-atom radial distribution functions (RDFs), and the model represents the three sheets of the clay layer as a single layer of coarse-grained sites.²⁸ Many of the coarse-grain potentials of ion-sheet interactions in this model were derived using potential of mean force (PMF) matching with fixed sheets, and the various coarse-grain potentials were converged sequentially.²⁸ When deriving the potential between free ions and the bulk sites on the clay sheet, Suter and co-workers used the layer-averaged z-coordinate. This method faithfully reproduced the behavior of ions in the fluidized system but would likely not be sufficiently accurate for reproducing the behavior of confined ions, which is strongly influenced by charge localization.

Simulations with both DFT and ClayFF indicate that ions in the anhydrous illite interlayer approach interlayer-facing oxygen atoms extremely closely, suggesting that compensatory compression of the neighboring sheets could play an important role in determining the energy landscape.²⁵ In addition, potassium and cesium ions at equilibrium within the interlayer are tightly confined within ditrigonal coordination cavities with oxygen, with interatomic distances of roughly 3.0–3.4 Å.^8,19 Due to this extreme confinement and the detailed structure of ion binding sites, representing an anhydrous clay layer by a single sheet of coarse-grained sites would likely eliminate important stiff stretching and bending degrees of freedom for determining the energetics of ion diffusion.

We present a coarse-grained model of an anhydrous illite clay system with different coarse-grain types for representing both the tetrahedral and octahedral sheets. This model runs approximately 70 times faster than the all-atom implementation in ClayFF within the LAMMPS⁴² simulation environment due to both reduced number of interaction sites as well as permitting much larger time steps of up to 10 fs for molecular dynamics (MD) simulations. While ClayFF is typically used to study only a few layers, our CG model can easily simulate dozens of clay layers. We show that studying large systems is necessary to completely eliminate finite size effects in the determination of converged diffusion barriers and provide evidence of a thermodynamic compensation for interstratification of potassium and cesium ion distributions in anhydrous illite clay interlayers. Our CG model is capable of investigating large systems (10⁻⁶ m) on simulation time scales of microseconds, and the model is available to others through the LAMMPS simulation package.

Classical molecular dynamics (MD) simulations of an atomistic 2-layer illite clay system under periodic boundary conditions were used as a reference for generating a coarse-grained (CG) model using iterative Boltzmann inversion. The all-atom MD simulations were performed using the ClayFF⁴³ force field within the LAMMPS⁴² simulation package. The ClayFF force field is a generalized, nonbonded model for hydrated clays, and consists mainly of Lennard-Jones and electrostatic interactions between atomic centers for bulk clay. It has been fitted to multiple multi-sheet, aluminosilicate clay types including kaolinite (Al₂Si₂O₅(OH)₄) and pyrophyllite (AlSi₂O₅(OH)) and has been shown to reproduce the swelling behavior of montmorillonite (Na₃(Si₃₁Al)(Al₁₄Mg₂)O₈₀(OH)₁₆) very accurately.⁴³ In addition, ClayFF has been used extensively to study the dynamics of ion adsorption in hydrated interlayers.^44–46 Because of its ability to model multilayered clay systems under a variety of physical conditions, ClayFF was chosen as the reference atomistic force field for our coarse-grained model.

All-atom molecular dynamics simulations were run at 300 K for 100 ps with a time step of 1 fs after an initial equilibration period of 120 ps. Coordinates of all atoms were sampled every 250 fs to build an ensemble for using the IBI algorithm. Each coarse-grain simulation was sampled in the NVT ensemble for 40 ps after an equilibration of 100 ps with a time step of 3 fs to compromise between fast turnaround time and sufficient sampling for IBI. To ensure faithful reproduction of the all-atom data, the coarse-grained systems used during the IBI procedure were mapped directly from the corresponding all-atom systems.

Climbing-Image Nudged Elastic Band (CINEB) calculations were performed using the “neb” command within LAMMPS^47–50 to obtain energy barriers for interlayer counterion (K⁺ and Cs⁺) migration, using 25 images integrated with a 5 fs time step. Before beginning each NEB calculation, both ion position and substitution sites were randomized. The first image for each NEB calculation was generated by equilibrating the randomized structure, and the final image had one ion from the first image displaced to an empty binding site. The remaining images were generated by linear interpolation of the ion position, so that the ion’s initial trajectory was a linear path between an occupied and unoccupied site. Oxygen atoms surrounding the initial and final ion binding sites were included in the reaction coordinate due to their displacement during ion diffusion. By including these atoms in the reaction coordinate, the distribution of energy barriers was reduced by around 30% without significantly altering the mean diffusion barrier value, indicating that this accurately captures important physical effects.

To determine system enthalpies, five-interlayer clay systems with different patterns of ion interstratification, periodic in the x- and y-directions with dimensions of 93 Å and 60 Å, respectively, were relaxed in the non-periodic z-direction over a period of 150 ps using a 3 fs time step. The z-direction was non-periodic, and was initialized at 50 Å. The simulations were integrated using the multilevel summation method (MSM) real-space electrostatics.^51,52 This electrostatic integration method computes short-range interactions exactly and computes long-range interactions by decomposing the potential into a sum of smooth potentials that are integrated with a series of progressively coarser meshes. MSM has a competitive level of accuracy as the particle-mesh Ewald (PME) technique for calculating long-range electrostatic interactions and unlike PME can be used for non-periodic systems. Using the “shrink-wrap” feature in LAMMPS, the z-dimension of the simulation box was allowed to dynamically change over the course of each simulation until convergence. This procedure resulted in the system reaching the interlayer spacing that minimized the system enthalpy.

The iterative Boltzmann inversion (IBI) algorithm for coarse-graining attempts to reproduce all-atom pair correlation functions by constructing a pairwise interaction potential, a consequence of Henderson’s uniqueness theorem.^53,54 CG simulations are run iteratively, and the new pair correlation function g_CG(r) is used to update the old interaction potential. Assuming only pairwise effects, the all-atom pair correlation function g_A(r) of gases can be roughly approximated:

where A is an arbitrary constant, β is the thermodynamic beta, and u(r) is the potential as a function of the radial separation. The iterative Boltzmann inversion algorithm updates the coarse-grain potential based on the all-atom and CG pair correlation functions:

Henderson’s uniqueness theorem is only precisely true for a homogenous fluid, but IBI is a robust method of coarse-graining that works for heterogeneous fluids and other phases as well. In particular, good performance of the IBI algorithm can be expected if run simultaneously for “orthogonal” degrees of freedom, such as the inter- and intra-sheet forces in our CG model of illite. The iterative Boltzmann inversion procedure was interfaced with LAMMPS by updating non-bonded and bonded coefficients as well as tabulated potentials after each iteration.

The CG model of anhydrous clay consists of 5 different coarse-grain centers. This mapping reduces the number of tracked centers in ClayFF by approximately 2:1 for the bulk clay. The five ClayCG centers are as follows: 1 CG type for the octahedral sheet (type Al) corresponding to the structural Al³⁺ cations; 2 CG types for the tetrahedral sheet, representing O²⁻ anions directly adjacent to the interlayer, but are separated into CG types corresponding to oxygen anions near sites with and without isomorphic substitutions of Al³⁺ for Si⁴⁺ (type Os and O), so that their interactions with other CG centers are computed separately in order to capture differences in the binding site characteristics. In addition, there are two types of ions (type K and Cs). Silicon and aluminum are not explicitly tracked in the tetrahedral sheet. See Fig. 1 for a representation of the coarse-grained model in Visual Molecular Dynamics (VMD) software.⁵⁵ 

The CG sites experience five types of forces based on the following interaction potential: harmonic bonds and angles between neighboring centers in the octahedral and tetrahedral sheets, Lennard-Jones interactions between CG centers, electrostatics between ion and substituted-oxygen CG centers, and finally tabulated forces [Eq. (3)],

The Lennard-Jones interactions are defined between each pair of CG centers, and are not based on mixing rules; this approach was chosen because the CG centers represent different numbers of atoms from the all-atom model, and therefore, the pairwise interactions are unlikely to be characterized by effective radii. Lennard-Jones interactions are slightly less computationally expensive than tabulated forces and are sufficient for the purposes of this study for characterizing the forces between sheets. Tabulated potentials are used between the ion and tetrahedral CG centers to more accurately capture the nature of the binding sites. Lennard-Jones interactions and electrostatics are excluded for 2nd, 3rd, and 4th neighbors based on the bond and angle topology. The interaction between ions is governed only by electrostatics, since the average distance between ions in the simulation (approximately 7 Å) is large enough that the Lennard-Jones interactions from the all-atom model will be negligible. After each round of IBI, the updated CG potential is fitted to either a harmonic or Lennard-Jones function for most of the degrees of freedom, and this fit is used in the next round of MD simulation. The CG parameters were found by running the IBI algorithm in pure K-illite and Cs-illite. Additional details can be found in the supplementary material.

In some cases, the same parameters were used to characterize multiple interactions due to the corresponding all-atom model centers having very similar pair correlation functions. For example, one single Lennard-Jones potential was used to govern the nonbonded interaction between the octahedral and tetrahedral CG sites. This approach further simplifies the model, although it does not affect its computational cost. To run our CG clay model, all that is needed is a properly configured input script and data file with the clay system coordinates, charges, and topology. No modification to LAMMPS software itself is needed to run the model, and the software used to run the IBI algorithm is entirely separate from the core codebase of any simulation software. The CG model is highly portable and should be able to run on any other molecular simulation package that allows for the implementation of tabulated potentials. The CG parameters are tabulated in the supplementary material, and the tabulated potentials are available in the LAMMPS format as supplementary material.

Figure 2 and Fig. S1 of the supplementary material show the distribution of bond distances and angles for the all-atom and CG models for the degrees of freedom fit with harmonic bonds, in which we observe overall excellent agreement. For the non-bonded degrees of freedom, Fig. 3 and Fig. S2 of the supplementary material present the radial distribution function (RDF) of the all-atom model in comparison to the most converged iteration for the CG model. Since the fit for the nonbonded degrees of freedom to the Lennard-Jones functions used only the first peak of each RDF, there is a relatively good agreement for all pair correlations with the exception of the O–O and O–Al CG types. However, even in these cases there is relatively good reproduction of the positions of secondary peaks in the RDFs, indicating that the geometry of the all-atom and CG systems is quite similar. The model is well converged based on the similarity of the oxygen-ion CG center RDF to the corresponding distribution in the all-atom model, since these degrees of freedom are the most important for fully characterizing the ion binding site.

In order to quantify the increase in speed for our coarse-grained model, an all-atom model and its corresponding coarse-grained model were run for 3 ns on 32 cores, and the real time needed to simulate each 100 fs was recorded. The all-atom system consisted of 2 clay layers with periodic boundary conditions, and a total of 9099 atoms and was simulated using the ClayFF force field. The coarse-grained model was simulated using the force field developed in this paper. The coarse-grained model ran roughly 6.9 times faster than the corresponding all-atom model with the same time step. The shortest vibrational period in the all-atom model is on the order of 10 fs due to the explicit modeling of hydrogen, and by contrast, the fastest vibrational mode in the coarse-grained model presented here is between oxygen-oxygen centers, which is on the order of 100 fs.⁵⁶ Because of this, the coarse-grained model is stable with a time step of up to 10 fs, while ClayFF must use an integration time step on the order of 1 fs. Thus our coarse-grained model represents a roughly 70-fold speedup compared to the corresponding all-atom force field, which is comparable to the speedup obtained by other groups using similar coarse-graining techniques.^40,41

As shown by Johnson et al.,⁵⁷ the Henderson uniqueness proof implies that there is always a representability problem as a general feature of a CG potential, i.e., a CG procedure cannot simultaneously resolve all the properties at a given state point. For example, reproducing the energetics of a system when the coarse-graining approach is based on reproducing structural or geometric features of the more complex reference system is not formally guaranteed.^20,58 However, the structural coarse-graining approach is likely to reproduce qualitative trends in properties such as energetic barriers for ion diffusion in anhydrous clays, since these barriers will be primarily determined by mechanical forces.

We sought to further validate our CG model by performing NEB calculations to determine diffusion barriers for different ions in the presence of the same or different ions in the clay interlayers. While the energy barriers derived from the coarse-grain model were consistently higher and had a broader distribution than the corresponding barriers in the all-atom model, Fig. 4 confirms that the CG model qualitatively reproduces the difference in diffusion energy barriers for K⁺ and Cs⁺ found previously in the all-atom clay model.²⁵ 

The mean energy barrier for the migration of potassium in pure K-illite was found to be 300 ± 94 kJ/mol [Fig. 4(a)] on average, compared to 226 ± 51 kJ/mol in the corresponding ClayFF model [Fig. 4(b)]. We also found for both K⁺ and Cs⁺, the energy barrier for diffusion was found to be much lower in systems with a higher fraction of Cs⁺ atoms in the interlayer. This result is in agreement with the trend derived from the all-atom force field, which has led to our determination of a mechanism for interlayer exchange,²⁵ whereby as more Cs⁺ enters an interlayer, both ion species become much more mobile, effectively increasing the rate at which the exchange front will advance into the interlayer.^7,8,28

Both the CG and all-atom energy barriers correspond to time scales that are inaccessible for direct observation of diffusion events in MD simulations,²⁵ but the ability of the CG model to approximate the energy barriers and their trends with respect to all-atom ClayFF is promising for using this coarse-grain model as a probe for changes in diffusion barriers to infer the kinetics of ion exchange. Properly modeling ion diffusion near and far from the exchange front necessarily requires a model capable of simulating a large, heterogeneous interlayer, as well as overcoming finite size effects by modeling many interlayers (i.e., greater than 2), which is extremely computationally expensive in all-atom ClayFF. Coarse-graining typically smoothens the energy landscape,⁵⁹ so it would be expected that the energy barriers in our CG model should be lower than the barriers in the all-atom model. However, the IBI methodology is relatively insensitive to the structure of the pairwise potential in undersampled regions, which are inherent in highly structured materials (see Fig. 3). The relatively large energy barriers in the ion–oxygen pair potentials (Fig. S4 of the supplementary material) are therefore likely an artifact of using the IBI procedure in a confined setting and are likely the source of the increased energy barriers.

For sufficiently small periodic systems, the finite size effect can dramatically impact the compressibility and would be expected to increase the calculated barrier to ion diffusion artificially.⁶⁰ Therefore, the magnitude of the NEB barriers was studied as a function of the number of simulated layers using periodic 2-interlayer, 4-interlayer, and 12-interlayer systems. Table I shows the average energy barrier for K⁺ diffusion in K-illite as a function of the number of clay layers, and the energy barrier distributions determined by NEB are presented in Fig. S3 of the supplementary material. Since there was essentially no observed change in the average energy barrier and variance between 4- and 12-interlayer systems, it was determined that simulating at least 4 interlayers would be sufficient to approximate an effectively infinite clay for the purposes of this study.

Mixing of unlike ions in layered silicate interlayers can adopt different structures depending on ion distributions in the final structure. When interlayer ions form random mixtures, the phase is homostructured, and when ions are separated into distinct phase-separated layers, the phase is interstratified [Figs. 5 and 6(a), respectively]. In the following section, we analyze the impact of the structure on the barriers to ion migration and consequently on the kinetics of ion exchange.

The results in Table II summarize the energy barriers to ion migration as a function of layer Cs⁺ content for homostructured K/Cs-illites (Fig. 5), in which each interlayer consisted of both Cs⁺ and K⁺ ions positioned randomly at counterion binding sites. As expected, both ions experience a statistically significant reduction in average energy barrier in Cs-illite compared to K-illite due to interlayer expansion, and there is a nearly linear trend in the change in barrier height as a function of composition, consistent with prior results using an all-atom force field.²⁵ The energy barrier distributions are quite broad, and there is significant overlap between the NEB barrier distributions in K-illite and Cs-illite. As a function of the change in equilibrium interlayer spacing, these results correspond to a decrease in the average barrier energy of approximately −71 (kJ/mol)/Å for Cs⁺ and −78 (kJ/mol)/Å for K⁺. This result is similar in magnitude to the change in energy barrier found in the all-atom ClayFF model of −92 (kJ/mol)/Å for K⁺.²⁵ The information in Table II is presented graphically in Fig. 8.

One advantage of our ClayCG model is the ability to model numerous clay structures with significantly reduced computational cost. Small clay systems were used to investigate the mixing enthalpy of interstratification due to ion exchange. These systems consisted of four periodic anhydrous interlayer regions between five illite clay layers stacked in a vertical configuration. The outer layers lack counterions on the exterior basal surfaces, which is necessary to allow convergence of the simulation cell size during the shrinkwrap procedure in LAMMPS. Although redistributing these ions in the interlayers is not physically realistic, it is not expected to significantly alter the equilibrium interlayer spacing or NEB energies presented in Table III, because basal spacing is controlled primarily by the counterion size. The excess mixing enthalpy (ΔH_mix) was computed as the deviation from ideal mixing between K- and Cs-illite,

where H_X is the computed minimum enthalpy of the clay system being simulated, f_i is the fraction of ion i in the clay system, and H_Cs and H_K are the enthalpies of five-layer clay containing only cesium and potassium in their interlayers, respectively (see Sec. II). For all of the five-layer systems investigated, each interlayer was occupied exclusively by one type of ion and did not have an exchange front. In the following tables, each five-layer system is abbreviated using the identity of the ions in its interlayers from the bottom to the top as a code. For example, the system corresponding to “Cs Cs K K” had two adjacent interlayers filled with cesium ions below two interlayers filled with potassium ions (see Fig. 6).

Table III summarizes the trends in the layer basal spacing, energy barriers to migration, and normal stress in the z-direction on both types of ions within different interstratification arrangements. For both K⁺ and Cs⁺ ions, the overall trend is towards increasing normal stress (and therefore cohesion) with decreasing interlayer spacing. In general, the thermodynamic compensation ΔH_mix is greatest in systems with more Cs-illite/K-illite interfaces. However, in the case of the Cs K Cs Cs particle, the expansion of the K-illite interlayer is not sufficient to compensate for compression of the Cs-illite interlayer. This result indicates that there is a penalty associated with disrupting contiguous Cs-illite particles.

The presence of Cs-illite decreases z-axial stress on K-illite layers, while the presence of K⁺ in the structure tends to increase z-axial stress on Cs⁺ interlayers. This finding suggests that the longstanding supposition² that exchanges on a layer increases cohesive energy in neighboring K⁺ interlayers is incorrect, at least when the adjacent layer also contains anhydrous counterions. Instead, we find that exchange on one layer, regardless of its proximity to the clay interior, alters the basal spacing and cohesive energies of ions in immediately adjacent layers. As shown in Fig. 7, the magnitude of this effect decreases rapidly as a function of distance from the exchanged layers and is nearly undetectable only a few layers away from the Cs-illite/K-illite interface.

For both types of ions, there is a clear trend of increasing energy barrier under increased confinement. In the case of only one Cs-illite layer present in bulk K-illite, it is clear from the interlayer spacing and normal stress that the cesium interlayer experiences maximal compression. To confirm this, NEB calculations were run on the displacement of a Cs⁺ ion in Cs-illite in the middle of a 12-layer K-illite particle. The barrier in this case was found to be 288 kJ/mol, very similar to the 286 kJ/mol barrier for Cs⁺ diffusion in the “K Cs K K” particle (Table III), supporting the conclusion that interlayer compression and expansion in interstratified particles is a localized effect. Similarly, a K-illite layer isolated in bulk cesium relaxes to a much larger spacing. To confirm that the results in Table III were not artifacts of small particle size, the average interlayer spacing in much larger particles (11 interlayers) of pure Cs-illite and K-illite were measured. These were found to be 1.053 ± 0.023 and 0.987 ± 0.022 nm, respectively, consistent with the results above.

The compensatory expansion and compression of neighboring layers may explain why the change in Cs⁺ and K⁺ diffusion barriers in interstratified particles is greater than the corresponding change in homostructured particles for a given interlayer spacing (Fig. 8). In homostructured clays, the local interlayer spacing near a Cs⁺ counterion is greater than the average interlayer spacing due to its larger atomic radius, resulting in a lower migration barrier for a given spacing. The opposite is true for K⁺ counterions, which experience greater local confinement in homostructured clays than would be expected from measuring the average interlayer spacing alone. By contrast, the compression and expansion of interlayers in interstratified particles results in a more uniform interlayer spacing and therefore a more dramatic change in the diffusion barrier. This effect may be an artifact of using ClayFF as the all-atom model, as ClayFF is known to slightly overestimate flexibility in the out-of-plane direction in large clay systems.⁶¹ 

From the trends in the energy barriers to ion diffusion in K-illite interlayers adjacent to Cs-illite interlayers (Table III), one would assume that exchange would be enhanced near cesium-dominated interlayers instead of being inhibited. That is, diffusion energy barriers in neighboring K-illite interlayers are lowered in the immediate vicinity of a Cs-illite interlayer. Experimental evidence for the Cs/K system is insufficient to confirm or refute this hypothesis, but there is some visual evidence that Cs exchanged layers occur in clumps in a K-phlogopite and are not randomly distributed, as expected from these results.⁹ 

Exchange of ions of different sizes leads to bending deformation of the layer structure locally, which may alter the coordination of K⁺ in the vicinity of the exchange front.⁶² In this case, the selectivity for and mobility of K⁺ is expected to vary with the sharpness and uniformity of the front. In order to capture the effects of the sharp exchange front in anhydrous illite clay systems,^7,8,25 a series of simulations were performed in which all ions on one half of the exchanging layer were assigned to be Cs⁺ and all ions on the other half were assigned to be K⁺. The exchange front was modeled in a periodic, four-layer clay system that featured one completely exchanged interlayer below the exchange front and two fully unexchanged interlayers above the exchange front. Ions far from the exchange front did not have significantly altered energy barriers to diffusion compared to barriers in an interstratified particle, but K⁺ ions characterized at the interface showed a slightly reduced average energy barrier compared to an interstratified clay with no exchange front barrier (273 kJ/mol vs. 278 kJ/mol). In comparison with the results of the energy barrier distribution presented in Table II, K⁺ ions near the exchange front experience approximately 20% additional barrier lowering (27 kJ/mol) with respect to the barrier in pure K-illite as compared to ions far from the front (22 kJ/mol), and ions more than approximately 2–3 nm from the exchange front are essentially unaffected by it. The exchange front was studied in an “intermediate” swelling state, adjacent to both Cs-illite and K-illite interlayers. The result indicates that the presence of an exchange front significantly alters potassium ion energy barriers even in the presence of an adjacent Cs-illite interlayer.

Amongst all of the non-periodic four interlayer systems studied, the system “K Cs K Cs” exhibited the greatest thermodynamic compensation for interstratification. Overall, the ΔH_mix trends indicate that there is a thermodynamic driver during the exchange process to form alternating K-illite and Cs-illite layers (i.e., ordered interstratified structures).^11,29,63 Exchange will tend to disrupt bulk K-illite, and thermodynamic feedback will favor exchange that leads to the formation of ordered interstratification instead of regions of bulk Cs-illite. However, the thermodynamic compensation for forming ordered interstratification is quite small (at most around 10 kJ/mol), indicating that it is unlikely to be the primary cause of experimentally observed differences in the exchange rate of adjacent interlayers.^7,8 Moreover, this compensation is significantly smaller than measured differences in hydration-free energies of aqueous K⁺ and Cs⁺,⁶⁴ which is expected to dominate the thermodynamics of ion exchange in hydrated clays and therefore near the frayed edge of illite clays.⁶⁵ It is expected that at least within a few nanometers of the clay edge, ordered interstratification is unlikely to be observed.

We have presented a coarse-grained model of anhydrous K- and Cs-illite which represents a 70-fold speedup over the corresponding all-atom ClayFF force field. Although using tabulated potentials between the oxygen coarse-grain centers may slightly improve the fidelity of modeling the ion binding sites, there is a non-negligible computational cost savings associated with using Lennard-Jones potentials as compared to splined tabulated potentials.⁶⁶ While other CG clay models have represented ion binding sites as single CG centers,^40,41 the model presented here is capable of capturing physical degrees of freedom important during ion diffusion in the confined interlayer.^8,19 By modeling all the atoms within and adjacent to the interlayer, we were able to accurately reproduce the structure of ion binding sites without significant computational overhead.

The reduction in particle density due to the coarse-grain procedure is especially helpful for nudged elastic band simulations, since the number of molecular dynamics steps necessary to reach convergence of the energy pathway is heavily dependent on the system size and the number of particles surrounding the transition path. We were able to converge NEB pathways on relatively large systems, which is promising for probing the energy barriers in physical scenarios that would necessitate investigating bulk effects. Our CG model qualitatively reproduced the ion diffusion energy trends as a function of interlayer separation found in ClayFF²⁵ and demonstrates interlayer expansion near sharp exchange fronts and near fully exchanged layers. Our model indicates a significant enthalpy of mixing associated with adjacent K- and Cs-illite interlayers in interstratified particles, which may impact exchange front propagation in adjacent interlayers. The coarse-graining strategy used here is expected to generalize well to other anhydrous and swelling clay types, and the intrasheet CG bond potentials will be directly portable to other clay systems.

See supplementary material for the parameters of the coarse-grained model and additional NEB results.

The authors acknowledge support by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.