Article pubs.acs.org/JPCA
Characterization of Adsorbed Alkali Metal Ions in 2:1 Type Clay Minerals from First-Principles Metadynamics Takashi Ikeda,*,† Shinichi Suzuki,‡,§ and Tsuyoshi Yaita‡,§ †
Condensed Matter Science Division and ‡Energy and Environment Materials Science Division, Quantum Beam Science Center, Japan Atomic Energy Agency (JAEA), 1-1-1 Kouto, Sayo-cho, Sayo-gun, Hyogo 679-5148, Japan § Fukushima Environmental Safety Center, Japan Atomic Energy Agency (JAEA), 6-6 Sakae-machi, Fukushima-shi, Fukushima 960-8031, Japan ABSTRACT: Adsorption states of alkali metal ions in three kinds of 2:1 type clay minerals are systematically investigated via first-principles-based metadynamics. Our reconstructed free energy surfaces in a two-dimensional space of coordination numbers specifically employed as collective variables for describing the interlayer cations show that an inner-sphere (IS) complex is preferentially formed for Cs+ in the 2:1 type trioctahedral clay minerals with saponite-like compositions, where lighter alkali metal ions show a tendency to form an outer-sphere one instead. The strong preference for an IS complex observed for Cs+ is found to result partially from the capability of recognizing selectively Cs+ ions at the basal O atoms with the Lewis basicity significantly enhanced by the isomorphic substitution in tetrahedral sheets.
■
INTRODUCTION Adsorption states of interlayer cations in swelling clay minerals such as the smectite and vermiculite groups are difficult to characterize quantitatively using experimental techniques particularly at the atomic level. Hence, the structure of water and the distribution of counterions in the interlayer of clay minerals have been addressed by using molecular simulations. Many simulation studies found in the literature suggest that adsorbed cations form either an inner-sphere (IS) or outersphere (OS) complex.1−10 In the former complex, interlayer cations are directly associated with isomorphic substitution (ISub) sites in tetrahedral (Th) sheets; thus the resulting cations are strongly bound to the clay surface. In contrast, in the latter one a cation is fully hydrated in the midplane of the interlayer, leading to the formation of electrical double layers in the interlayer in a complicated fashion. For instance, the free energy profile for the adsorption of a sodium ion onto the surface of a smectite clay evaluated by using constrained firstprinciples molecular dynamics (FPMD)5 shows that the lowest free-energy region corresponds to a Na+ bound to only one basal O (Ob) atom, accompanied by the shift of the position of a Na+ ion by 1.5 Å from the midplane of the interlayer toward the Th sheet including the ISub. In contrast, it is suggested from FPMD that Na+ ions in Na vermiculite are located at the midplane of the interlayer, though the interlayer contents appear to be disordered.7 Thus, the adsorption state of interlayer cations is likely to vary depending on details of clay minerals including the location of ISub sites, the amount of layer charge, and the activity of interlayer water. However, to the best of our knowledge, fundamental physicochemistry responsible for the behavior of adsorbed ions and coexisting © XXXX American Chemical Society
water in swelling clay minerals is not yet fully understood, though clay swelling plays a key role in various applications such as oil and gas production.11 This work is motivated by recent experimental findings related to the Cs+ adsorption in clay minerals. Motokawa et al.12 suggest from small-angle X-ray scattering (SAXS) that the abruption of clay layers is induced by the collective adsorption of Cs+ ions in the interlayer. Kogure et al.13 directly observed the fixation of Cs+ in the interlayer of vermiculite clay using high-resolution transmission electron microscopy (HRTEM). Moreover, radiocesium isotopes 134Cs and 137Cs released in the Fukushima Daiichi nuclear disaster occurring in 2011 were observed to aggregate in particular clay minerals of weathered biotite constituting soils around Fukushima.14 However, detailed mechanisms for aggregating Cs+ ions in the interlayer of particular clay minerals are not yet understood to date. In this computational study, we describe the detailed structure of adsorbed alkali metal ions in 2:1 type clay minerals by reconstructing the free energy surface (FES). Although the structures of our clay systems including the ISub in either Th or octahedral (Oh) sheets only are oversimplified compared to the actual ones, the present systematic study reveals the underlying fundamental sciences in varied adsorption states of a series of alkali metal ions in the interlayer of our model clays, thus providing a clue to understanding the cesium aggregation in particular clay minerals. Received: June 21, 2015 Revised: July 6, 2015
A
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
■
COMPUTATIONAL METHODS Setup for Metadynamics. The FES in the space spanned by a set of collective variables (CVs) selected specifically for describing the adsorbed interlayer cations was obtained by performing metadynamics simulations.15,16 In this work, we employed two coordination numbers (CNs) as CVs, namely, the metal (M)−Ob CN nMOb and the M−water O (Ow) CN nMOw defined as16,17 nMO{s} =
1 NM
NM
∑∑ j ∈ {s}
i
assigned to our CVs and spring constant for restraining of k = 0.24 au ensure the adiabatic separation of our CVs from the other degrees of freedom. To estimate the FES with high resolution, small Gaussians with a height of H < 0.31 kcal mol−1 and width of W = 0.1 were added to the history-dependent potential included in the formulation of metadynamics according to a distance criterion of 3/2 W as recommended by Ensing et al.21 The electronic structure was computed by using revPBE GGA functional22,23 complemented by the first-principles van der Waals corrections (denoted as vdW-WF) based on the use of maximally localized Wannier centers of the type reported in refs 24 and 25. Our specific implementation in the CPMD code26 is documented in the code’s manual and will be reported in detail elsewhere. The valence−core interaction was described by a norm-conserving Troullier−Martins27 pseudopotential (PP) for H, O, Na, Mg, Al, Si, and K, and Goedecker28,29 PP for Li, Rb, and Cs generated based on the full relativistic electronic structure calculations of a single atom. Valence orbitals were expanded in a plane-wave (PW) basis set with an energy cutoff of 80 Ry and periodic boundary conditions on the simulation cell were applied. The Brillouin zone was sampled at the Γ point only. Our PP-based PW-DFT approach was validated by examining the ion−oxygen distances d and related hydration energies ΔE for M+H2O complexes, as demonstrated in refs 30 and 18. Our estimated d and ΔE differ at most by 0.06 Å and 2.3 kcal mol−1, respectively, from the corresponding estimates made with a standard quantum chemistry method at the MP2 level of theory.30−32 Preparation of Our Clay Systems. Our systems were prepared by performing constant pressure Born−Oppenheimer MD simulations with variable cell shape33 for the 2:1 type clay minerals at P = 1 atm and T = 300 K. The equations of motion were integrated with a time step of 20.671 au (0.5 fs). The electronic structure was computed by using the same combination of revPBE GGA functional and the vdW-WF as in metadynamics simulations. We used an energy cutoff of 140 Ry to ensure the convergence of the computed stress tensor within the accuracy of 10 kbar. Figure 1 shows snapshots of our three kinds of clays taken from the trajectories generated in our preparation runs. Clay minerals have crystal structures consisting of stacked phyllosilicates composed of Th and Oh sheets. Each clay layer of the 2:1 type clay minerals is composed of two Th sheets consisting dominantly of SiO4 units fused to an Oh sheet as schematically represented in Figure 2. Note that all our clays contain 24 H2O molecules in the interlayer, which are sufficient for creating a water bilayer, as shown in Figure 1. As shown in Figure 1a, our saponite (denoted as Sap), belonging to the tri-Oh subgroup, includes the ISub of Al3+ for
1 − (rM iOj /rc)n 1 − (rM iOj /rc)n + m
(1)
where NM is the number of metal atoms, rMiOj is the interatomic distance between the ith M and jth O atoms, and the symbol {s} represents a subset of O atoms in our clay systems represented below. To use CNs as CVs of metadynamics, we should carefully determine the value of a cutoff distance rc so that the computed CN value using eq 1 is as close as possible to the (integer) number of oxygen atoms included in the coordination shell of a given metal ion. Our selected values of rc for defining the two CNs nMOb and nMOw are tabulated in Table 1. Note that our selected values of rc for MOw pairs are Table 1. Values of a Cutoff Distance rc (Å) Selected for Defining Two CNs nMOb and nMOw pair
Li
Na
K
Rb
Cs
MOb MOw
3.30 2.90
3.90 3.40
4.40 3.70
4.45 3.80
4.50 3.90
0.4−0.7 Å smaller than those for MOb ones. This is due to our previous observation that the boundary between the first and second solvation shells of heavy alkali metal ions in bulk water becomes ambiguous.18 In this particular application the two exponents n and m controlling the decay of the rational function in eq 1 around rc were selected as 10 and 12, respectively. In this study metadynamics simulations were performed within the Car−Parrinello framework19 based on the density functional theory20 (DFT) under constant NVT conditions at T = 300 K for the three kinds of 2:1 type clay minerals containing a series of alkali metal ions in the interlayer by employing the averaged lattice parameters generated in our preparation runs described below. The temperature was controlled by rescaling of ionic velocities. The equations of motion were integrated with a time step of 3 au (0.0726 fs) and a fictitious electron mass of 340 au, ensuring reliable static and dynamical properties. A large fictitious mass of M = 20 amu
Figure 1. Snapshots of our (a) Cs-Sap, (b) Cs-Bei, and (c) Cs-Mnt taken from preparation runs performed under constant NPT conditions of 1 atm and 300 K. The instantaneous simulation cell is shown as black solid lines. Atom colors are white for H, red for O, cyan for Mg, pink for Al, yellow for Si, and silver for Cs. B
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
lengths of 13.8−15.9 Å are in reasonable agreement with the available experimental ones. Therefore, we expect that our clay systems reproduce characteristic features of the 2:1 type clay minerals even though their structures are oversimplified compared to the actual ones.
■
RESULTS AND DISCUSSION Reconstructed FESs. Figure 3 shows contour plots of our reconstructed FESs as a function of the two CNs nMOb and nMOw for our Sap clays containing alkali metal ions in their interlayer. In this simulation each 3 Ob atoms next to an Al3+ ion in each Th sheet were used for computing nMOb. Consequently, the computed value of nMOb ranges up to 3.0 from its definition and our selected rc values given in eq 1 and Table 1, respectively. Typical configurations corresponding to the free-energy minima labeled as A, B, ..., H in Figure 3 are shown in Figure 4. For Li+ the minimum A located at around (0.0, 4.1) in Figure 3a corresponds to the OS complex, in which 4 H2O molecules participate (Figure 4a), as observed in bulk water30 (and in our Bei and Mnt clays). Our FES shows that the metastable IS complex corresponding to the local minimum at (1.5, 3.2) in Figure 3a is ∼5 kcal mol−1 higher than the most stable OS tetra-aqua complex. For Na+ the region of the minimum B in Figure 3b, corresponding to an OS complex as represented in Figure 4b, is broader than that of the minimum A for Li+, indicating that an OS complex of Na+ is more flexible than that of Li+. Interestingly, our FES for K+ shows three minima at around (0.2, 6.2), (1.1, 5.4), and (1.8, 5.5) (Figure 3c), where each K+ ion tends to be bound to at least one Ob atom in Th sheets, as shown in Figure 4c,d. Although K+ ions have a tendency to form an IS complex in the interlayer rather than an OS one, the presence of the three shallow minima indicates that their affinity with ISub sites in Th sheets is not so high. Our FES for Rb+ (Figure 3d) also shows three minima at around (0.5, 5.6), (1.1, 5.1), and (1.9, 5.0), where an IS complex is dominantly formed (Figure 4e,f), as observed for K+. However, the relative stability of the three minima for Rb+ is shown to be reversed compared to that for K+, indicating that the affinity of Rb+ with ISub sites in Th sheets is much higher than that of K+. For Cs+ two minima G and H in Figure 3e are found to be separated from each other by a rather high activation barrier of ∼5 kcal mol−1. The former is a local minimum, corresponding to Cs+ ions being not bound to any Ob atoms next to Al3+ ions (Figure 4g), whereas the latter is a global minimum, corresponding to the IS complex associated with 2 Ob atoms next to an ISub site of Al3+ (Figure 4h). Therefore, our results suggest that Cs+ ions strongly prefer to form the particular IS complex associated with ISub sites in our Sap. Figure 5 shows the reconstructed FESs for our Bei and Mnt clays including two Na+, K+, and Cs+ ions in their interlayer. Here, each 3 Ob atoms next to an Al3+ (a selected Si4+) in each Th sheet were used for computing nMOb for our Bei (Mnt) clays. As shown in Figure 5a, the FES for Na+ shows a shallow minimum at around (0.8, 5.1) along with a global one at (0.1, 5.3) in our Bei. The position of the latter minimum for Na+ is observed to shift from (0.1, 5.3) (Figure 5a) to (1.1, 5.1) (Figure 5d) by changing ISub sites from Th sheets as in our Bei to Oh ones as in our Mnt. On the contrary, for K+ the region of the minimum extending in the range 0 Mnt among our clays, in good accord with that of the observed affinity of Cs+ with Ob atoms. Note that for Na+ the order of the affinity with Ob atoms appears to be inverted, in line with our observations that interlayer H2O molecules rather than Na+ ions tend to work as a Lewis acid for interacting with Ob atoms.
■
CONCLUSIONS In conclusion, we have successfully obtained the FESs for adsorption states of alkali metal ions in the interlayer of the three representative 2:1 type clay minerals via FPMD-based metadynamics. Cs+ ions are found to be selectively recognized among alkali metal ions at the Ob atoms in Th sheets of the 2:1 type tri-Oh clays with Sap-like compositions, where the Lewis basicity of the Ob atoms is significantly enhanced by the ISub of Al3+ for their neighboring Si4+ in Th sheets. Because for lighter alkali metal ions interlayer H2O molecules rather than the ions themselves tend to work as a Lewis acid for interacting with such Ob atoms, lighter alkali metal ions show a tendency to form an OS complex instead of an IS one preferred by heavier alkali metal ions. Furthermore, our close examination of electronic structures indicates that the Cs−O interactions in the IS complex preferentially formed for Cs+ ions in Sap-like clays include small covalency, thus leading to significantly high affinity of Cs+ with those Ob atoms. In contrast, in our Bei and Mnt clays the Lewis basicity of Ob atoms is significantly reduced particularly in our Mnt compared to that in our Sap clays. Therefore, the characteristic adsorption state for light and heavy alkali metal ions observed in our Sap tends to disappear in our Bei and Mnt clays except for Li+. Hence, the 2:1 type clay minerals suitably constituted are expected to work as a superior extractor of Cs+, which is highly likely to be pertinent to the Cs+ aggregation.
■
AUTHOR INFORMATION
Corresponding Author
*T. Ikeda. Phone: +81-791-58-2622. Fax: +81-791-58-0311. Email: [email protected]. Notes
The authors declare no competing financial interest.
■
Figure 7. Computed LDOSs at Cs, Ob, and Al sites along with energyresolved overlap population QCsOb(E) for our Sap clay including a Cs+ ion bound to 2 Ob atoms next to an Al3+ in a Th sheet.
ACKNOWLEDGMENTS This work was performed under the JAEA project, “Cs sorption-desorption mechanism on clay minerals”, based on the special account for Fukushima environment recovery from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The computational work reported in this paper was performed using the supercomputing facilities in JAEA and Information Technology Center of The University of Tokyo.
Figure 7 shows the LDOSs and energy-resolved overlap population QCsOb(E) for our Sap clay containing a Cs+ ion bound to 2 Ob atoms next to an Al3+ ion in a Th sheet. Compared to the LDOS of Ob 2p shown in Figure 6a, the corresponding LDOS shown in Figure 7 is observed to be significantly broadened, indicating a non-negligible contribution of chemical bonding between Ob 2p and some orbitals of a Cs+ ion upon the cesium adsorption. The QCsOb(E) for the resulting Cs−Ob bond displayed in the bottom panel of Figure 7 indicates that the contribution of the bonding orbitals between Cs+ 5sp and Ob 2p is almost completely canceled by that of the counterpart antibonding orbitals widely distributed from −3.9 to +2.0 eV, though the contribution of bonding orbitals at −4.3 eV in which Cs+ 5p, Ob 2p, and Al 3s participate seems to be large. Thus, the net contribution to the chemical bonding between Cs+ and Ob atoms results mainly from the electron donation from Ob 2p to Cs+ 5d empty orbitals. Indeed, the LDOS of Cs+ 5d shows some small features in the occupied region, as displayed in the inset of Figure 7. Therefore, the present analysis suggests that a unique electronic feature of
■
REFERENCES
(1) Sposito, G.; Skipper, N. T.; Sutton, R.; Park, S.-h.; Soper, A. K.; Greathouse, J. A. Surface Geochemistry of the Clay Minerals. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 3358−3364. (2) Marry, V.; Turq, P. Microscopic Simulations of Interlayer Structure and Dynamics in Bihydrated Heteroionic Montmorillonites. J. Phys. Chem. B 2003, 107, 1832−1839. (3) Boek, E. S.; Sprik, M. Ab Initio Molecular Dynamics Study of the Hydration of a Sodium Smectite Clay. J. Phys. Chem. B 2003, 107, 3251−3256. (4) Whitley, H. D.; Smith, D. E. Free Energy, Energy, and Entropy of Swelling in Cs-, Na-, and Sr-Montmorillonite Clays. J. Chem. Phys. 2004, 120, 5387−5395. (5) Suter, J. L.; Boek, E. S.; Sprik, M. Adsorption of a Sodium Ion on a Smectite Clay from Constrained Ab Initio Molecular Dynamics Simulations. J. Phys. Chem. C 2008, 112, 18832−18839. F
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A (6) Mignon, P.; Ugliengo, P.; Sodupe, M.; Hernandez, E. R. Ab Initio Molecular Dynamics Study of the Hydration of Li+, Na+ and K+ in a Montmorillonite Model. Influence of Isomorphic Substitution. Phys. Chem. Chem. Phys. 2010, 12, 688−697. (7) Demontis, P.; Masia, M.; Suffritti, G. B. Water Nanoconfined in Clays: The Structure of Na Vermiculite Revisited by Ab Initio Simulations. J. Phys. Chem. C 2013, 117, 15583−15592. (8) Demontis, P.; Masia, M.; Suffritti, G. B. Peculiar Structure of Water in Slightly Superhydrated Vermiculite Clay Studied by CarParrinello Molecular Dynamics Simulations. J. Phys. Chem. C 2014, 118, 7923−7931. (9) 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. (10) Sun, L.; Tanskanen, J. T.; Hirvi, J. T.; Kasa, S.; Schatz, T.; Pakkanen, T. A. Molecular Dynamics Study of Montmorillonite Crystalline Swelling: Roles of Interlayer Cation Species and Water Content. Chem. Phys. 2015, 455, 23−31. (11) Valaškova, M., Martynkova, G. S., Eds. Clay Minerals in Nature Their Characterization, Modification and Application; InTech: Rijeka, Croatia, 2012. (12) Motokawa, R.; Endo, H.; Yokoyama, S.; Nishitsuji, S.; Kobayashi, T.; Suzuki, S.; Yaita, T. Collective Structural Changes in Vermiculite Clay Suspensions Induced by Cesium Ions. Sci. Rep. 2014, 4, 6585. (13) Kogure, T.; Morimoto, K.; Tamura, K.; Sato, H.; Yamagishi, A. XRD and HRTEM Evidence for Fixation of Cesium Ions in Vermiculite Clay. Chem. Lett. 2012, 41, 380−382. (14) Mukai, H.; Hatta, T.; Kitazawa, H.; Yamada, H.; Yaita, T.; Kogure, T. Speciation of Radioactive Soil Particles in the Fukushima Contaminated Area by IP Autoradiography and Microanalyses. Environ. Sci. Technol. 2014, 48, 13053−13059. (15) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (16) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302. (17) Stirling, A.; Iannuzzi, M.; Laio, A.; Parrinello, M. Azulene-toNaphthalene Rearrangement: The Car-Parrinello Metadynamics Method Explores Various Reaction Mechanisms. ChemPhysChem 2004, 5, 1558−1568. (18) Ikeda, T.; Boero, M. Communication: Hydration Structure and Polarization of Heavy Alkali Ions: A First Principles Molecular Dynamics Study of Rb+ and Cs+. J. Chem. Phys. 2012, 137, 041101. (19) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (20) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (21) Ensing, B.; Laio, A.; Parrinello, M.; Klein, M. L. A Recipe for the Computation of the Free Energy Barrier and the Lowest Free Energy Path of Concerted Reactions. J. Phys. Chem. B 2005, 109, 6676−6687. (22) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (23) Zhang, Y.; Yang, W. Comment on ”Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1998, 80, 890. (24) Silvestrelli, P. Van der Waals Interactions in DFT Made Easy by Wannier Functions. Phys. Rev. Lett. 2008, 100, 053002. (25) Silvestrelli, P. L. Van der Waals Interactions in Density Functional Theory Using Wannier Functions. J. Phys. Chem. A 2009, 113, 5224−5234. (26) CPMD, version 3.17.1; The CPMD consortium, www.cpmd.org, 2013. (27) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for PlaneWave Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006.
(28) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (29) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 3641−3662. (30) Ikeda, T.; Boero, M.; Terakura, K. Hydration of Alkali Ions from First Principles Molecular Dynamics Revisited. J. Chem. Phys. 2007, 126, 034501. (31) Kołaski, M.; Lee, H. M.; Choi, Y. C.; Kim, K. S.; Tarakeshwar, P.; Miller, D. J.; Lisy, J. M. Structures, Energetics, and Spectra of AquaCesium (I) Complexes: An Ab Initio and Experimental Study. J. Chem. Phys. 2007, 126, 074302. (32) San-Román, M.; Hernández-Cobos, J.; Saint-Martin, H.; OrtegaBlake, I. A Theoretical Study of the Hydration of Rb+ by Monte Carlo Simulations with Refined Ab Initio-Based Model Potentials. Theor. Chem. Acc. 2010, 126, 197−211. (33) Focher, P.; Chiarotti, G. L.; Bernasconi, M.; Tosatti, E.; Parrinello, M. Structural Phase Transformations via First-Principles Simulation. Europhys. Lett. 1994, 26, 345.
G
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Characterization of Adsorbed Alkali Metal Ions in 2:1 Type Clay Minerals from First-Principles Metadynamics Takashi Ikeda,*,† Shinichi Suzuki,‡,§ and Tsuyoshi Yaita‡,§ †
Condensed Matter Science Division and ‡Energy and Environment Materials Science Division, Quantum Beam Science Center, Japan Atomic Energy Agency (JAEA), 1-1-1 Kouto, Sayo-cho, Sayo-gun, Hyogo 679-5148, Japan § Fukushima Environmental Safety Center, Japan Atomic Energy Agency (JAEA), 6-6 Sakae-machi, Fukushima-shi, Fukushima 960-8031, Japan ABSTRACT: Adsorption states of alkali metal ions in three kinds of 2:1 type clay minerals are systematically investigated via first-principles-based metadynamics. Our reconstructed free energy surfaces in a two-dimensional space of coordination numbers specifically employed as collective variables for describing the interlayer cations show that an inner-sphere (IS) complex is preferentially formed for Cs+ in the 2:1 type trioctahedral clay minerals with saponite-like compositions, where lighter alkali metal ions show a tendency to form an outer-sphere one instead. The strong preference for an IS complex observed for Cs+ is found to result partially from the capability of recognizing selectively Cs+ ions at the basal O atoms with the Lewis basicity significantly enhanced by the isomorphic substitution in tetrahedral sheets.
■
INTRODUCTION Adsorption states of interlayer cations in swelling clay minerals such as the smectite and vermiculite groups are difficult to characterize quantitatively using experimental techniques particularly at the atomic level. Hence, the structure of water and the distribution of counterions in the interlayer of clay minerals have been addressed by using molecular simulations. Many simulation studies found in the literature suggest that adsorbed cations form either an inner-sphere (IS) or outersphere (OS) complex.1−10 In the former complex, interlayer cations are directly associated with isomorphic substitution (ISub) sites in tetrahedral (Th) sheets; thus the resulting cations are strongly bound to the clay surface. In contrast, in the latter one a cation is fully hydrated in the midplane of the interlayer, leading to the formation of electrical double layers in the interlayer in a complicated fashion. For instance, the free energy profile for the adsorption of a sodium ion onto the surface of a smectite clay evaluated by using constrained firstprinciples molecular dynamics (FPMD)5 shows that the lowest free-energy region corresponds to a Na+ bound to only one basal O (Ob) atom, accompanied by the shift of the position of a Na+ ion by 1.5 Å from the midplane of the interlayer toward the Th sheet including the ISub. In contrast, it is suggested from FPMD that Na+ ions in Na vermiculite are located at the midplane of the interlayer, though the interlayer contents appear to be disordered.7 Thus, the adsorption state of interlayer cations is likely to vary depending on details of clay minerals including the location of ISub sites, the amount of layer charge, and the activity of interlayer water. However, to the best of our knowledge, fundamental physicochemistry responsible for the behavior of adsorbed ions and coexisting © XXXX American Chemical Society
water in swelling clay minerals is not yet fully understood, though clay swelling plays a key role in various applications such as oil and gas production.11 This work is motivated by recent experimental findings related to the Cs+ adsorption in clay minerals. Motokawa et al.12 suggest from small-angle X-ray scattering (SAXS) that the abruption of clay layers is induced by the collective adsorption of Cs+ ions in the interlayer. Kogure et al.13 directly observed the fixation of Cs+ in the interlayer of vermiculite clay using high-resolution transmission electron microscopy (HRTEM). Moreover, radiocesium isotopes 134Cs and 137Cs released in the Fukushima Daiichi nuclear disaster occurring in 2011 were observed to aggregate in particular clay minerals of weathered biotite constituting soils around Fukushima.14 However, detailed mechanisms for aggregating Cs+ ions in the interlayer of particular clay minerals are not yet understood to date. In this computational study, we describe the detailed structure of adsorbed alkali metal ions in 2:1 type clay minerals by reconstructing the free energy surface (FES). Although the structures of our clay systems including the ISub in either Th or octahedral (Oh) sheets only are oversimplified compared to the actual ones, the present systematic study reveals the underlying fundamental sciences in varied adsorption states of a series of alkali metal ions in the interlayer of our model clays, thus providing a clue to understanding the cesium aggregation in particular clay minerals. Received: June 21, 2015 Revised: July 6, 2015
A
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
■
COMPUTATIONAL METHODS Setup for Metadynamics. The FES in the space spanned by a set of collective variables (CVs) selected specifically for describing the adsorbed interlayer cations was obtained by performing metadynamics simulations.15,16 In this work, we employed two coordination numbers (CNs) as CVs, namely, the metal (M)−Ob CN nMOb and the M−water O (Ow) CN nMOw defined as16,17 nMO{s} =
1 NM
NM
∑∑ j ∈ {s}
i
assigned to our CVs and spring constant for restraining of k = 0.24 au ensure the adiabatic separation of our CVs from the other degrees of freedom. To estimate the FES with high resolution, small Gaussians with a height of H < 0.31 kcal mol−1 and width of W = 0.1 were added to the history-dependent potential included in the formulation of metadynamics according to a distance criterion of 3/2 W as recommended by Ensing et al.21 The electronic structure was computed by using revPBE GGA functional22,23 complemented by the first-principles van der Waals corrections (denoted as vdW-WF) based on the use of maximally localized Wannier centers of the type reported in refs 24 and 25. Our specific implementation in the CPMD code26 is documented in the code’s manual and will be reported in detail elsewhere. The valence−core interaction was described by a norm-conserving Troullier−Martins27 pseudopotential (PP) for H, O, Na, Mg, Al, Si, and K, and Goedecker28,29 PP for Li, Rb, and Cs generated based on the full relativistic electronic structure calculations of a single atom. Valence orbitals were expanded in a plane-wave (PW) basis set with an energy cutoff of 80 Ry and periodic boundary conditions on the simulation cell were applied. The Brillouin zone was sampled at the Γ point only. Our PP-based PW-DFT approach was validated by examining the ion−oxygen distances d and related hydration energies ΔE for M+H2O complexes, as demonstrated in refs 30 and 18. Our estimated d and ΔE differ at most by 0.06 Å and 2.3 kcal mol−1, respectively, from the corresponding estimates made with a standard quantum chemistry method at the MP2 level of theory.30−32 Preparation of Our Clay Systems. Our systems were prepared by performing constant pressure Born−Oppenheimer MD simulations with variable cell shape33 for the 2:1 type clay minerals at P = 1 atm and T = 300 K. The equations of motion were integrated with a time step of 20.671 au (0.5 fs). The electronic structure was computed by using the same combination of revPBE GGA functional and the vdW-WF as in metadynamics simulations. We used an energy cutoff of 140 Ry to ensure the convergence of the computed stress tensor within the accuracy of 10 kbar. Figure 1 shows snapshots of our three kinds of clays taken from the trajectories generated in our preparation runs. Clay minerals have crystal structures consisting of stacked phyllosilicates composed of Th and Oh sheets. Each clay layer of the 2:1 type clay minerals is composed of two Th sheets consisting dominantly of SiO4 units fused to an Oh sheet as schematically represented in Figure 2. Note that all our clays contain 24 H2O molecules in the interlayer, which are sufficient for creating a water bilayer, as shown in Figure 1. As shown in Figure 1a, our saponite (denoted as Sap), belonging to the tri-Oh subgroup, includes the ISub of Al3+ for
1 − (rM iOj /rc)n 1 − (rM iOj /rc)n + m
(1)
where NM is the number of metal atoms, rMiOj is the interatomic distance between the ith M and jth O atoms, and the symbol {s} represents a subset of O atoms in our clay systems represented below. To use CNs as CVs of metadynamics, we should carefully determine the value of a cutoff distance rc so that the computed CN value using eq 1 is as close as possible to the (integer) number of oxygen atoms included in the coordination shell of a given metal ion. Our selected values of rc for defining the two CNs nMOb and nMOw are tabulated in Table 1. Note that our selected values of rc for MOw pairs are Table 1. Values of a Cutoff Distance rc (Å) Selected for Defining Two CNs nMOb and nMOw pair
Li
Na
K
Rb
Cs
MOb MOw
3.30 2.90
3.90 3.40
4.40 3.70
4.45 3.80
4.50 3.90
0.4−0.7 Å smaller than those for MOb ones. This is due to our previous observation that the boundary between the first and second solvation shells of heavy alkali metal ions in bulk water becomes ambiguous.18 In this particular application the two exponents n and m controlling the decay of the rational function in eq 1 around rc were selected as 10 and 12, respectively. In this study metadynamics simulations were performed within the Car−Parrinello framework19 based on the density functional theory20 (DFT) under constant NVT conditions at T = 300 K for the three kinds of 2:1 type clay minerals containing a series of alkali metal ions in the interlayer by employing the averaged lattice parameters generated in our preparation runs described below. The temperature was controlled by rescaling of ionic velocities. The equations of motion were integrated with a time step of 3 au (0.0726 fs) and a fictitious electron mass of 340 au, ensuring reliable static and dynamical properties. A large fictitious mass of M = 20 amu
Figure 1. Snapshots of our (a) Cs-Sap, (b) Cs-Bei, and (c) Cs-Mnt taken from preparation runs performed under constant NPT conditions of 1 atm and 300 K. The instantaneous simulation cell is shown as black solid lines. Atom colors are white for H, red for O, cyan for Mg, pink for Al, yellow for Si, and silver for Cs. B
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
lengths of 13.8−15.9 Å are in reasonable agreement with the available experimental ones. Therefore, we expect that our clay systems reproduce characteristic features of the 2:1 type clay minerals even though their structures are oversimplified compared to the actual ones.
■
RESULTS AND DISCUSSION Reconstructed FESs. Figure 3 shows contour plots of our reconstructed FESs as a function of the two CNs nMOb and nMOw for our Sap clays containing alkali metal ions in their interlayer. In this simulation each 3 Ob atoms next to an Al3+ ion in each Th sheet were used for computing nMOb. Consequently, the computed value of nMOb ranges up to 3.0 from its definition and our selected rc values given in eq 1 and Table 1, respectively. Typical configurations corresponding to the free-energy minima labeled as A, B, ..., H in Figure 3 are shown in Figure 4. For Li+ the minimum A located at around (0.0, 4.1) in Figure 3a corresponds to the OS complex, in which 4 H2O molecules participate (Figure 4a), as observed in bulk water30 (and in our Bei and Mnt clays). Our FES shows that the metastable IS complex corresponding to the local minimum at (1.5, 3.2) in Figure 3a is ∼5 kcal mol−1 higher than the most stable OS tetra-aqua complex. For Na+ the region of the minimum B in Figure 3b, corresponding to an OS complex as represented in Figure 4b, is broader than that of the minimum A for Li+, indicating that an OS complex of Na+ is more flexible than that of Li+. Interestingly, our FES for K+ shows three minima at around (0.2, 6.2), (1.1, 5.4), and (1.8, 5.5) (Figure 3c), where each K+ ion tends to be bound to at least one Ob atom in Th sheets, as shown in Figure 4c,d. Although K+ ions have a tendency to form an IS complex in the interlayer rather than an OS one, the presence of the three shallow minima indicates that their affinity with ISub sites in Th sheets is not so high. Our FES for Rb+ (Figure 3d) also shows three minima at around (0.5, 5.6), (1.1, 5.1), and (1.9, 5.0), where an IS complex is dominantly formed (Figure 4e,f), as observed for K+. However, the relative stability of the three minima for Rb+ is shown to be reversed compared to that for K+, indicating that the affinity of Rb+ with ISub sites in Th sheets is much higher than that of K+. For Cs+ two minima G and H in Figure 3e are found to be separated from each other by a rather high activation barrier of ∼5 kcal mol−1. The former is a local minimum, corresponding to Cs+ ions being not bound to any Ob atoms next to Al3+ ions (Figure 4g), whereas the latter is a global minimum, corresponding to the IS complex associated with 2 Ob atoms next to an ISub site of Al3+ (Figure 4h). Therefore, our results suggest that Cs+ ions strongly prefer to form the particular IS complex associated with ISub sites in our Sap. Figure 5 shows the reconstructed FESs for our Bei and Mnt clays including two Na+, K+, and Cs+ ions in their interlayer. Here, each 3 Ob atoms next to an Al3+ (a selected Si4+) in each Th sheet were used for computing nMOb for our Bei (Mnt) clays. As shown in Figure 5a, the FES for Na+ shows a shallow minimum at around (0.8, 5.1) along with a global one at (0.1, 5.3) in our Bei. The position of the latter minimum for Na+ is observed to shift from (0.1, 5.3) (Figure 5a) to (1.1, 5.1) (Figure 5d) by changing ISub sites from Th sheets as in our Bei to Oh ones as in our Mnt. On the contrary, for K+ the region of the minimum extending in the range 0 Mnt among our clays, in good accord with that of the observed affinity of Cs+ with Ob atoms. Note that for Na+ the order of the affinity with Ob atoms appears to be inverted, in line with our observations that interlayer H2O molecules rather than Na+ ions tend to work as a Lewis acid for interacting with Ob atoms.
■
CONCLUSIONS In conclusion, we have successfully obtained the FESs for adsorption states of alkali metal ions in the interlayer of the three representative 2:1 type clay minerals via FPMD-based metadynamics. Cs+ ions are found to be selectively recognized among alkali metal ions at the Ob atoms in Th sheets of the 2:1 type tri-Oh clays with Sap-like compositions, where the Lewis basicity of the Ob atoms is significantly enhanced by the ISub of Al3+ for their neighboring Si4+ in Th sheets. Because for lighter alkali metal ions interlayer H2O molecules rather than the ions themselves tend to work as a Lewis acid for interacting with such Ob atoms, lighter alkali metal ions show a tendency to form an OS complex instead of an IS one preferred by heavier alkali metal ions. Furthermore, our close examination of electronic structures indicates that the Cs−O interactions in the IS complex preferentially formed for Cs+ ions in Sap-like clays include small covalency, thus leading to significantly high affinity of Cs+ with those Ob atoms. In contrast, in our Bei and Mnt clays the Lewis basicity of Ob atoms is significantly reduced particularly in our Mnt compared to that in our Sap clays. Therefore, the characteristic adsorption state for light and heavy alkali metal ions observed in our Sap tends to disappear in our Bei and Mnt clays except for Li+. Hence, the 2:1 type clay minerals suitably constituted are expected to work as a superior extractor of Cs+, which is highly likely to be pertinent to the Cs+ aggregation.
■
AUTHOR INFORMATION
Corresponding Author
*T. Ikeda. Phone: +81-791-58-2622. Fax: +81-791-58-0311. Email: [email protected]. Notes
The authors declare no competing financial interest.
■
Figure 7. Computed LDOSs at Cs, Ob, and Al sites along with energyresolved overlap population QCsOb(E) for our Sap clay including a Cs+ ion bound to 2 Ob atoms next to an Al3+ in a Th sheet.
ACKNOWLEDGMENTS This work was performed under the JAEA project, “Cs sorption-desorption mechanism on clay minerals”, based on the special account for Fukushima environment recovery from the Ministry of Education, Culture, Sports, Science and Technology of Japan. The computational work reported in this paper was performed using the supercomputing facilities in JAEA and Information Technology Center of The University of Tokyo.
Figure 7 shows the LDOSs and energy-resolved overlap population QCsOb(E) for our Sap clay containing a Cs+ ion bound to 2 Ob atoms next to an Al3+ ion in a Th sheet. Compared to the LDOS of Ob 2p shown in Figure 6a, the corresponding LDOS shown in Figure 7 is observed to be significantly broadened, indicating a non-negligible contribution of chemical bonding between Ob 2p and some orbitals of a Cs+ ion upon the cesium adsorption. The QCsOb(E) for the resulting Cs−Ob bond displayed in the bottom panel of Figure 7 indicates that the contribution of the bonding orbitals between Cs+ 5sp and Ob 2p is almost completely canceled by that of the counterpart antibonding orbitals widely distributed from −3.9 to +2.0 eV, though the contribution of bonding orbitals at −4.3 eV in which Cs+ 5p, Ob 2p, and Al 3s participate seems to be large. Thus, the net contribution to the chemical bonding between Cs+ and Ob atoms results mainly from the electron donation from Ob 2p to Cs+ 5d empty orbitals. Indeed, the LDOS of Cs+ 5d shows some small features in the occupied region, as displayed in the inset of Figure 7. Therefore, the present analysis suggests that a unique electronic feature of
■
REFERENCES
(1) Sposito, G.; Skipper, N. T.; Sutton, R.; Park, S.-h.; Soper, A. K.; Greathouse, J. A. Surface Geochemistry of the Clay Minerals. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 3358−3364. (2) Marry, V.; Turq, P. Microscopic Simulations of Interlayer Structure and Dynamics in Bihydrated Heteroionic Montmorillonites. J. Phys. Chem. B 2003, 107, 1832−1839. (3) Boek, E. S.; Sprik, M. Ab Initio Molecular Dynamics Study of the Hydration of a Sodium Smectite Clay. J. Phys. Chem. B 2003, 107, 3251−3256. (4) Whitley, H. D.; Smith, D. E. Free Energy, Energy, and Entropy of Swelling in Cs-, Na-, and Sr-Montmorillonite Clays. J. Chem. Phys. 2004, 120, 5387−5395. (5) Suter, J. L.; Boek, E. S.; Sprik, M. Adsorption of a Sodium Ion on a Smectite Clay from Constrained Ab Initio Molecular Dynamics Simulations. J. Phys. Chem. C 2008, 112, 18832−18839. F
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A (6) Mignon, P.; Ugliengo, P.; Sodupe, M.; Hernandez, E. R. Ab Initio Molecular Dynamics Study of the Hydration of Li+, Na+ and K+ in a Montmorillonite Model. Influence of Isomorphic Substitution. Phys. Chem. Chem. Phys. 2010, 12, 688−697. (7) Demontis, P.; Masia, M.; Suffritti, G. B. Water Nanoconfined in Clays: The Structure of Na Vermiculite Revisited by Ab Initio Simulations. J. Phys. Chem. C 2013, 117, 15583−15592. (8) Demontis, P.; Masia, M.; Suffritti, G. B. Peculiar Structure of Water in Slightly Superhydrated Vermiculite Clay Studied by CarParrinello Molecular Dynamics Simulations. J. Phys. Chem. C 2014, 118, 7923−7931. (9) 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. (10) Sun, L.; Tanskanen, J. T.; Hirvi, J. T.; Kasa, S.; Schatz, T.; Pakkanen, T. A. Molecular Dynamics Study of Montmorillonite Crystalline Swelling: Roles of Interlayer Cation Species and Water Content. Chem. Phys. 2015, 455, 23−31. (11) Valaškova, M., Martynkova, G. S., Eds. Clay Minerals in Nature Their Characterization, Modification and Application; InTech: Rijeka, Croatia, 2012. (12) Motokawa, R.; Endo, H.; Yokoyama, S.; Nishitsuji, S.; Kobayashi, T.; Suzuki, S.; Yaita, T. Collective Structural Changes in Vermiculite Clay Suspensions Induced by Cesium Ions. Sci. Rep. 2014, 4, 6585. (13) Kogure, T.; Morimoto, K.; Tamura, K.; Sato, H.; Yamagishi, A. XRD and HRTEM Evidence for Fixation of Cesium Ions in Vermiculite Clay. Chem. Lett. 2012, 41, 380−382. (14) Mukai, H.; Hatta, T.; Kitazawa, H.; Yamada, H.; Yaita, T.; Kogure, T. Speciation of Radioactive Soil Particles in the Fukushima Contaminated Area by IP Autoradiography and Microanalyses. Environ. Sci. Technol. 2014, 48, 13053−13059. (15) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (16) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302. (17) Stirling, A.; Iannuzzi, M.; Laio, A.; Parrinello, M. Azulene-toNaphthalene Rearrangement: The Car-Parrinello Metadynamics Method Explores Various Reaction Mechanisms. ChemPhysChem 2004, 5, 1558−1568. (18) Ikeda, T.; Boero, M. Communication: Hydration Structure and Polarization of Heavy Alkali Ions: A First Principles Molecular Dynamics Study of Rb+ and Cs+. J. Chem. Phys. 2012, 137, 041101. (19) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (20) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (21) Ensing, B.; Laio, A.; Parrinello, M.; Klein, M. L. A Recipe for the Computation of the Free Energy Barrier and the Lowest Free Energy Path of Concerted Reactions. J. Phys. Chem. B 2005, 109, 6676−6687. (22) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (23) Zhang, Y.; Yang, W. Comment on ”Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1998, 80, 890. (24) Silvestrelli, P. Van der Waals Interactions in DFT Made Easy by Wannier Functions. Phys. Rev. Lett. 2008, 100, 053002. (25) Silvestrelli, P. L. Van der Waals Interactions in Density Functional Theory Using Wannier Functions. J. Phys. Chem. A 2009, 113, 5224−5234. (26) CPMD, version 3.17.1; The CPMD consortium, www.cpmd.org, 2013. (27) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for PlaneWave Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006.
(28) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (29) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 3641−3662. (30) Ikeda, T.; Boero, M.; Terakura, K. Hydration of Alkali Ions from First Principles Molecular Dynamics Revisited. J. Chem. Phys. 2007, 126, 034501. (31) Kołaski, M.; Lee, H. M.; Choi, Y. C.; Kim, K. S.; Tarakeshwar, P.; Miller, D. J.; Lisy, J. M. Structures, Energetics, and Spectra of AquaCesium (I) Complexes: An Ab Initio and Experimental Study. J. Chem. Phys. 2007, 126, 074302. (32) San-Román, M.; Hernández-Cobos, J.; Saint-Martin, H.; OrtegaBlake, I. A Theoretical Study of the Hydration of Rb+ by Monte Carlo Simulations with Refined Ab Initio-Based Model Potentials. Theor. Chem. Acc. 2010, 126, 197−211. (33) Focher, P.; Chiarotti, G. L.; Bernasconi, M.; Tosatti, E.; Parrinello, M. Structural Phase Transformations via First-Principles Simulation. Europhys. Lett. 1994, 26, 345.
G
DOI: 10.1021/acs.jpca.5b05934 J. Phys. Chem. A XXXX, XXX, XXX−XXX
