18832
J. Phys. Chem. C 2008, 112, 18832–18839
Adsorption of a Sodium Ion on a Smectite Clay from Constrained Ab Initio Molecular Dynamics Simulations J. L. Suter* Department of Chemistry, UniVersity of Cambridge, Cambridge, CB2 1EW, United Kingdom and Schlumberger Cambridge Research, High Cross, Madingley Road, Cambridge, CB3 0EL, United Kingdom
E. S. Boek Schlumberger Cambridge Research, High Cross, Madingley Road, Cambridge, CB3 0EL, United Kingdom
M. Sprik Department of Chemistry, UniVersity of Cambridge, Cambridge, CB2 1EW, United Kingdom ReceiVed: July 27, 2007; ReVised Manuscript ReceiVed: June 29, 2008
We have performed ab intitio molecular dynamics (MD) simulations to study the adsorption of a sodium ion onto the surface of a smectite clay. In a previous study using Car-Parrinello MD simulations [Boek and Sprik, J. Phys. Chem. B 2003, 107, 3251], spontaneous adsorption was not observed. To study this process in more detail, we have used constraint MD simulations to sample regions of the phase space inaccessible to standard MD methods. We have chosen a reaction coordinate and constrained the value over a certain range to force adsorption. Using this method, we have calculated potentials of mean force for the process. We also have obtained, for the first time, the free-energy differences between a fully solvated “outer sphere” complex, and the sodium that is bound to the surface (an “inner sphere” complex). These complexes have been identified in previous experiments and simulations using classical potentials. We found that the lowest free-energy region exists at a distance of 6.1 Å away from the center of the clay layer, between the “inner-sphere” and “outer-sphere” complexes. In this state, the sodium is bound to only one surface oxygen atom. The stabilization of this region is greater when the sodium adsorbs onto a Si5Al ring. 1. Introduction The swelling of clays is a process that plays a role in many industrial processes, such as petroleum engineering and catalysis.4,5 Generally, the interactions between alkali-metal ions and mineral surfaces such as clays and carbonates are of fundamental interest for drilling fluids,6 as well as for enhanced oil recovery processes.7 Research has shown that 2:1 clays such as smectites and vermiculites consist of negatively charged aluminosilicate layers kept together by charge-balancing interlayer cations to give a stacked crystalline structure.6 Their most-characteristic property is the ability to adsorb water between the layers, hydrating the interlayer cations and resulting in strong repulsive forces and clay expansion. Additional water may then occupy the remaining space between the layers. To fully understand this swelling behavior, knowledge of the hydration of the counterions and the clay surface is required. Montmorillonites are commonly used as models for this effect, both in experiments and simulations.8-12 Sodium intercalates are the archetypal clay swelling materials.13 They are comprised of negatively charged pyrophyllite-like sheets, which are held together by chargebalancing sodium cations.6 The interlayer swelling of smectite minerals such as montmorillonite is believed to proceed stepwise through the formation of one-, two-, and three-layer hydrates.11,14,15 A two-layer hydrate, for example, can be pictured as two water * Author to whom correspondence should be addressed. Present address: Centre for Computational Science, Department of Chemistry, University College London, 20, Gordon Street, WC1H 0AJ, U.K. Fax: 442076791501. E-mail:
[email protected].
monolayers between the clay sheets.16 However, this process is dependent on the nature of the counterion; although Namontmorillonites are easily hydrated, Cs-montmorillonite is stopped at the monolayer, because of the low hydration capacity of cesium.17 The strong interaction of potassium (K+) ions with the clay surface in potassium smectite clays also results in the incomplete hydration of the counterion. This fact accounts for the use of potassium as a clay-swelling inhibitor.10 Quantitative details of the hydration behavior of montmorillonites are difficult to access experimentally, especially at the molecular level.18 The small crystallite size and the difficulty of producing wellordered aggregates hamper the study of the water structure in the interlayer region, and many counterions have small scattering intensities in neutron and X-ray diffraction. Investigations utilizing isotopic substitution have thrown some light on the structure of interlayer water, but so far, it has only been possible to analyze the coordination environment of counterions with large neutron scattering lengths, such as nickel (Ni+).16 Isotopic substitution (H/D and 6Li/natLi) of Li-montmorillonite hydrates have been undertaken, but it was not possible to determine whether the Li+ is partially or fully hydrated.19 Therefore, molecular simulation is an ideal way of gaining insight at the microscopic level.20-22 Simulation of the molecular structure of water and the counterion distribution of Na-, K-, and Li-montmorillonite clays has been addressed in previous studies using molecular dynamics (MD)11,12,23,24 and Monte Carlo (MC) techniques.10,14 Recently, electronic structure computer simulations of clays have appeared in the literature.22 Ab initio molecular dynamics (AIMD) simulations have shown
10.1021/jp075946a CCC: $40.75 2008 American Chemical Society Published on Web 11/07/2008
Adsorption of Na+ on Smectite Clay much promise, because this methodology removes the need for tuning pair-potentials to reproduce the characteristics of the clay system, such as the equilibrium swelling. Atomic forces are calculated by a continuously updated electronic structure calculation. The electronic structure calculation is based on the density functional theory (DFT) for periodic extended systems, as originally developed for application to solids. Structural optimization using DFT has been shown to be able to reproduce even subtle features of the pyrophyllite structure.26 Because of the coupling to MD, the finite temperature ionic dynamics is also described by this technique, thereby opening up the modeling of liquids to electronic structure methods. This allows us to account for the presence of the aqueous solvent, similar to classical simulation. This methodology for studying Na-montmorillonite clay using AIMD has been outlined in a previous study by Boek and Sprik, which will be called Paper I in the remainder of this study.1 These results confirmed that the water shows a strong preference to form a hydrogen-bonded network between the solvent molecules. The hydrogen bonds to the aluminosilicate surface were weak and short-lived. Unlike the ab initio MD study of muscovite mica, the solvent was shown to display liquid-like relaxation dynamics and the confinement prevents freezing, even in quantities as small as a double layer, albeit on the very limited time scale of the simulation.27-29 However, the rigor of AIMD comes at a considerable computational cost, with the accessible time scale being of the order of 1-10 ps, rather than the 1-10 ns expected of classical MD simulation. To circumvent this, in Paper I,1 a sample equilibrated snapshot from a classical MC simulation was adapted to fit the smaller periodicity of the AIMD simulation cell to form suitable initial conditions. A production run of 1 ps was performed, from which it was possible to observe the reorganization of the interlayer previously described. However, on the time scale of the simulation, the Na+ ion only explored the area around the middle of the interlayer spacing, near its location in the initial structure. From previous classical MD and MC simulations, we expect a two-layer hydrate with a Na+ counterion to possess two species, as defined by Sposito:30 an “inner-sphere complex”, where the sodium is strongly bound to the surface of the clay, associated with sites of substitution in the tetrahedral layer, and an “outer-sphere complex”, where the sodium is fully hydrated by water molecules in the midplane of the interlayer. This latter complex is associated with substitutions in the octahedral layer. In the study of Boek et al.,10 an “inner sphere” complex was observed with the sodium at a distance of approximately z ) 5.1 Å from the center of the clay layer (i.e., the octahedral ions), and it is this definition that we use in the remainder of the study. The absence of the “innersphere” complex in Paper I1 may be due to the short time scale of simulation that is possible with AIMD, especially as one would expect its prevalence because the unit-cell structure only contains isomorphic substitutions in the tetrahedral layer. Therefore, we require a method to allow the Na+ ion to sample the entire interlayer spacing, but on a time scale that is possible with AIMD. This is the focus of the current paper, which is a continuation of Paper I.1 To allow the Na+ ion to explore the interlayer spacing, we have decided to use Potential of Mean Force methods,2,3,31-33 which renders the study of Na+ ion adsorption possible. Using this method, we calculate potentials of mean force and the free-energy differences between the “outer-sphere” and “inner-sphere” complexes. This method involves a series of sequential simulations, each starting from the end of the
J. Phys. Chem. C, Vol. 112, No. 48, 2008 18833 previous simulation, at fixed constraint values. We define the constraint to be the distance in the z-direction between the Na+ ion and the center of the clay sheet. Although the length of each individual simulation within this series is short, using a large number of constraint values with only a small change in the constraint results in only a small perturbation to the system. The short simulation length possible for each constraint value means that the Na+ ion, although not constrained in the xyplane, may not explore all regions of the surface on which to adsorb. In an effort to counteract this, we also start the Na+ ion above a region of the clay sheet that does not contain an isomorphic substitution, and force adsorption onto this part of the clay sheet. This study is a first step to understanding the attraction of the clay surface to oppositely charged species. The classical MD/MC picture indicates that the surface is highly attractive and there is a fine balance between adsorption and hydration, leading to both “inner-sphere” and “outer-sphere” complexes. However, in Paper I,1 no “inner-sphere” complexes were formed and, indeed, hydrogen bonding to surface from interlayer water molecules was less than expected. Using constrained AIMD simulations, we attempt to describe the balance between adsorption and hydration for the Na+ ion. 2. Computational Methods 2.1. Clay Simulation Model. The simulation cell is identical to the cell dimensions in Paper I.1 The large size of unit cells of complex materials such as clay systems requires drastic simplifications to allow treatment by AIMD. The basis is the model Wyoming montmorillonite,18 which has been studied in previous classical Monte Carlo and MD simulations.10 This is a 2:1 layered dioctahedral aluminosilicate that belongs to the smectite family of clay materials and has the following unitcell formula:
Na0.75[Si7.75Al0.25][Al3.5Mg0.5]O20(OH)4 · nH2O
(1)
The simplifications required to allow structure 1 to be studied via AIMD involves using only two unit cells and one Na+ ion. This requires only one substitution in the clay framework. The substitution site we choose is an Si atom in the tetrahedral layer, which is replaced by an Al atom. Sixteen water molecules were added to this system, which is sufficient for the creation of a water double layer. These simplifications have been justified previously in Paper I.1 The rationale for preferring a tetrahedral substitution rather than an octahedral substitution is that it is presumed that the negative charge associated with the substitution is distributed only over the three surface O atoms bound to the defect Al atom. This should be more localized and less affected by periodic boundary conditions than in the octahedral case.30 Also, inner-sphere complexes, where the ion is bound to the surface of the clay, are normally assumed to be associated with tetrahedral substitution sites. Technically, the clay in the simulation cell is beidellite, because the only substitutions are in the tetrahedral layer. Both beidellite and montmorillonite are members of the smectite family. Therefore, the simulation box has the following formula:
Na1.00[Si15Al1][Al8]O40(OH)8 · 16H2O The dimensions of this cell are 10.56 Å × 9.14 Å, and with the water double layer, the z cell dimension is 15.01 Å. Periodic boundary conditions are applied in all three directions. In Figure 1, we present two simulations snapshots: one from the simulation, showing the Na+ ion in the middle of the interlayer (left;
18834 J. Phys. Chem. C, Vol. 112, No. 48, 2008
Figure 1. Simulation snapshot of a clay hydrate simulation. The simulation cell contains contains a total of 64 O, 9 Al, 15 Si, and 40 H atoms and 1 Na+ ion. This figure shows the simulation cell in the xz-plane. The left-hand panel shows the Na+ ion in the midplane of the interlayer at this constraint value, in an “outer-sphere complex”, fully hydrated by the interlayer water molecules. The right-hand panel shows the Na+ ion adsorbed on the surface in an “inner-sphere” complex, which is defined as being a distance of ∼5.1 Å away from the center of the clay layer (i.e., the Al3+ ions). It is bound to the surface of the clay near the site of aluminum substitution. Periodic boundary conditions are applied at the boundaries of the cell, indicated by the solid black lines.
an outer-sphere complex) and one absorbed on the surface (right; an inner-sphere complex). 2.2. Methodology. The AIMD scheme used is the CarParrinello method.25 This combines a quantum mechanical treatment of the electrons using the plane-wave pseudo-potentialbased density functional theory (DFT)34 method with a classical treatment of the nuclei, which are moved via an MD simulation. This scheme has been described in Paper I1 for Na-montmorillonite. The exchange functional of Becke35 and the correlation functional according to Lee, Yang, and Parr36 (known as the BLYP functional) were used. The Kohn-Sham orbitals are expanded in plane waves up to an energy cutoff of 70 Rydberg. The pseudo-potentials used are medium soft normconserving pseudo-potentials, according to the Troullier-Martins37 (TM) scheme for O, Si, H and the Bachelet, Hamann, and Schlu¨ter pseudo-potential for Al.38 Two different pseudo-potentials have been used for sodium. In the first case, the TM scheme was used, which includes the 2s and 2p states of Na as valence states, to properly address nonlinear exchange and correlation effects. In the second instance, a Bachelet, Hamann, and Schlu¨ter38,39-type pseudopotential is utilized, with only the 3s state considered as valence electrons. Details of the pseudo-potentials are given in Paper I1 and in an investigation of the properties of aqueous Na+ ions.40 Initial calculations were performed using the 2s,2p semicore potential for sodium. Although the much-smoother 3s,3p BHS potential previously, in combination with the BLYP functional, has given accurate structures and dynamics in aqueous solution, it was felt that when the sodium is bound to the surface of the clay, the semicore potential would be more reliable. Several test runs were performed with both potentials yielding little difference in the potentials of mean force, as shown in the Appendix. Therefore, for further calculations that involved the sodium adsorbed onto a Si6 ring on the clay surface, we have used the computationally less-demanding 3s BHS potential. The simulations were performed using a fictitious electron mass of 1000 au (1.66 × 10-24 kg) and a time step of 5 au (0.12 fs). All calculations were conducted under microcanonical conditions (NVT) with the lattice dimensions fixed during the simulations. The temperature was set at an initial value of 300 K, and velocity scaling ensured that the temperature was maintained within 100 K of this value in the equilibrium stages.
Suter et al. To calculate the potentials of mean force, we conducted a series of simulations with the z-coordinate of the Na+ ion fixed at various constraint values. The clay surface is located in the xy plane. The initial constraint value used in this study is the end-point of the simulation of Paper I, which had a production run of 1 ps.1 These simulations were performed sequentially; the end point of the previous constraint value is used as the starting point for the next constraint value with the Na+ ion moved manually to the next constraint value. The movement of the Na+ ion between simulations was never greater than 0.5 Å. Each run at a fixed constraint value lasted at least 1.5 ps. Equilibrium was assumed to be reached by 0.5 ps or until the electronic energy had reached a constant level, and data were collected from this point. We believe the short simulation times are valid, as the perturbation to the system between constraint values is small. One picosecond should be sufficient to accumulate enough data of interest, such as the forces required for the potential of mean force and structural details such as the radial distribution functions and atomic density profiles. On the time scale of our simulation, we see reorientation of the water molecules to accommodate the new constraint value, as shown by changes in the radial distribution functions. However, we do not calculate any time-dependent quantities, because of the short simulation timespan. To reduce any statistical error and allow greater sampling of the xy plane, one would need to run longer simulations. By gradually moving the Na+ ion onto the clay surface, the sodium was adsorbed onto a Si5Al ring. We also wished to calculate adsorption onto a Si6 ring of the clay surface. To achieve this, we chose a selected snapshot from the unconstrained simulation of Paper I where the sodium resided above a Si6 ring and manually moved the Na+ atom ∼1 Å closer to the surface. The wave function and geometry were optimized, and this served as a starting point for the Si6 adsorption constrained MD simulations. From the end of this simulation, the Na+ ion was then manually moved another 1.75 Å closer to the Si6 surface, at a z-coordinate that was consistent with an “inner-sphere” complex from previous classical MD simulations.10 Again, wave function and geometry optimizations were performed and this was used as a starting point for MD simulations. The remaining constraint values were constructed in a stepwise function from the equilibrium geometry of the nearest constraint value. In addition to the constrained MD, we have also performed two unconstrained MD simulations. First, we continued the simulation of Paper I with the Na+ ion in the midplane of the interlayer, up to a simulation length of 7.5 ps. A further unconstrained simulation used the inner-sphere complex, with a constraint value of 5.1 Å as an initial starting structure, but with the constraint removed. This simulation lasted 2 ps. Calculations were performed using the CPMD package.41,42 In total, ∼25 000 CPU h were used, with each simulation approaching 1000 CPU h. 3. Results and Discussion 3.1. Thermodynamics of Adsorption. The spontaneous adsorption of sodium onto the montmorillonite clay surface is not observable in a free AIMD run, because the reaction is an activated process. Therefore, we will study the reaction using constrained MD. The free-energy change for this process can be extracted by thermodynamic integration, provided that the control parameter is a reasonable approximation to the reaction coordinate. We have chosen the difference in the z-coordinates of the Na+ ion and the center of mass of the octahedral Al atoms
Adsorption of Na+ on Smectite Clay
J. Phys. Chem. C, Vol. 112, No. 48, 2008 18835
Figure 2. (a) Constraint force and (b) the corresponding change in free energy for the Na+ ion along the reaction coordinate. The dashed line is the force for adsorption onto a Si6 ring on the clay surface. The solid line is adsorption onto a Si5Al ring.
as the reaction coordinate, which is referred to as zNa-Al in the remainder of this study. The Na atom is allowed to move in the xy plane (i.e., move over the clay surface) only by constraining the difference in z-coordinates. At each point along the reaction coordinate, the ensembleaveraged force due to the constraint is evaluated. (This force is calculated as the Lagrange undetermined multiplier 〈f〉, which is determined by solving the equation of motion of the system with the constraint.) For the simple distance constraint that we are using in this study, the free-energy change is given by
dF ) - 〈f 〉 zNa-Al dzNa-Al
(2)
Therefore, the calculation requires a series of constrained MD runs at different values of the reaction coordinate. The freeenergy increments are obtained by integrating the average force along the reaction coordinate. As described in section 2.2, the constraint was varied between 9.9 Å and 4.3 Å in discrete steps. This samples the entire interlayer space. The first and last points are almost within the van der Waals radius of the clay surface O atoms, which lie at z ) 3.4 and 11.6 Å. The aluminum substitution in the tetrahedral layer lies at z ) 2.7 Å. The ensemble-averaged forces and free-energy curve are presented in Figure 2. A negative force, in this instance, indicates that the Na+ ion is being repelled away from the surface that contains the aluminum substitution toward the midplane of the interlayer spacing. We have conducted two simulations for adsorption onto different parts of the clay surface: the solid line in Figure 2 represents the Na+ ion lying above the Si5Al octahedral ring of the clay surface and the dashed curve indicates the Na+ adsorbing onto a Si6 ring of the clay surface. 3.2. Adsorption onto Si5Al Ring. At small constraint values (
