Water Nanoconfined in Clays: The Structure of Na Vermiculite

The results of Car–Parrinello molecular dynamics simulations of vermiculite clay at the hydration level found in the natural mineral (4 H2O per Na+)...
1 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCC

Water Nanoconfined in Clays: The Structure of Na Vermiculite Revisited by Ab Initio Simulations Pierfranco Demontis,† Marco Masia,†,§ and Giuseppe B. Suffritti*,† †

Dipartimento di Chimica e Farmacia, Università di Sassari and Consorzio Interuniversitario Nazionale per la Scienza e Tecnologia dei Materiali (INSTM), Unitá di ricerca di Sassari, Via Vienna, 2, I-07100 Sassari, Italy § Istituto Officina dei Materiali del CNR, UOS SLACS, Via Vienna 2, 07100 Sassari, Italy S Supporting Information *

ABSTRACT: The results of Car−Parrinello molecular dynamics simulations of vermiculite clay at the hydration level found in the natural mineral (4 H2O per Na+) are reported. The structure and the vibrational spectrum of the aluminosilicate layer is well reproduced by the simulations. As for the structure of the Na+ cations and of water molecules adsorbed in the interlayer, the ones proposed in X-ray diffraction experimental papers are not in full agreement with the results of neutron diffraction experiments found in literature, yielding only the density profile perpendicular to the aluminosilicate layers, but including hydrogens. Our calculations result in a structure of the interlayer content, which is in agreement with the neutron-diffraction density profile, and yet compatible with X-ray diffraction data. Some features of the adsorbed water resulting from the simulation are an example of incomplete cation hydration and can be of general interest.

1. INTRODUCTION Clay minerals are layer type aluminosilicates, which are ubiquitous on our planet in geologic deposits, terrestrial weathering environments, and marine sediments.1−5 The behavior of water confined in clays is an important clue for understanding many various phenomena in geophysics, catalysis, and environmental science fields. Clays are made of aluminosilicates layers (or sheets) of nanometric thickness, in which sheets of TO4 tetrahedra (usually T = Si, Al) are fused to sheets of MO4 octahedra (M = Al, Mg, Fe, Ti or other metallic cations with various electric charge). The layers may contain one tetrahedral and one octahedral sheet (1:1 layer type), or one octahedral sheet sandwiched between two tetrahedral sheets (2:1 layer type). Depending on the degree of isomorphic replacement of tetrahedral Si by Al, on the type and amount of octahedral cations, and on the number of hydroxyls that can form in the octahedra vertices, a negative charge is usually created in the sheets. It is balanced by interlayer cations (the most common are Na+, K+, and Ca2+). A variable amount of water molecules can be intercalated between clay layers to create an interlayer ionic solution that causes swelling phenomena related to electrical double layer properties and makes the interlayer cation exchange possible. The behavior of water depends on layer charge and on the type of interlayer cations. A widespread 2:1 layer type class of clays are smectites, whose layers have an average negative charge of about −0.5 e per unit cell (u.c.).1,4−6 Vermiculite too is a clay mineral belonging to the 2:1 layer type with negative structural charge −2 e per u.c., which is created by isomorphic Al/Si replacement. The interlayer cations are usually Na+ or Ca2+, and are adsorbed on or near © 2013 American Chemical Society

the basal plane of a tetrahedral sheet. Natural Na vermiculite contains about four water molecules per Na+ cation, but water can be removed by heating or adsorbed from a sufficiently humid environment. Its name comes from Latin vermiculare, to breed worms, for the manner in which it exfoliates when heated. It finds many different practical uses, mainly as insulating or fire protecting material, in packaging, and as substrate or additive in agricultural applications. In addition, the behavior of water in vermiculite clay of different temperatures and hydration levels has been investigated by several experimental techniques as a typical example of nanoconfined water.8−14 Therefore, a deep understanding of the vermiculite structure and in particular of the adsorbed water is of general physicochemical interest. The X-ray diffraction structure of Na vermiculite was reported by Slade et al.15 The positions of the atoms in the aluminosilicate sheet were given with good accuracy. Na+ cations were found in two partially occupied positions (m1 and m2), and the location of the water molecules, which must be coordinated to the cations and thus are subject to severe geometrical constraints, were not completely resolved. More recently, Beyer and Graf von Reichenbach,16 on the basis of the data reported in ref 15 and in other previous experimental studies, performed an extended revision of the vermiculite structure by applying geometrical constraints for atomic distances and bond angles in the data refinement. The structure of the interlayer hydrated cations was proposed as a chain-like arrangement of Na(H2O)6 octahedra, sharing two oxygen Received: January 23, 2013 Revised: July 1, 2013 Published: July 8, 2013 15583

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

were treated as valence states. The time step was set to 0.097 fs (4 au) and the fictitious mass for the orbital was chosen to be 400 amu (3.644 × 10−28 kg). This is the typical value used for simulations of bulk water with hydrogens.45 To check the validity of these choices, we performed test calculations for Na+(H2O)n clusters (n = 1 to 6) using the CPMD with the above basis set and parameters and ab initio DFT calculations with the def2-QZVP basis set. The PBE density functional was used in both cases. It can be seen in Figure S1 (Supporting Information) that the agreement is good. The simulated systems were first equilibrated at 330 K for 5 ps, using a massive Nosé-Hoover chain (NHC) thermostat;46−48 massive thermostatting consists of applying a NHC to each degree of freedom of the system, thus yielding to energy equipartition in a short time. It was then cooled at 290 K and the production simulations were run in the microcanonical ensemble. In the Supporting Information it is also shown that the total energy is well conserved and that the CPMD energy surface is parallel to the Born−Oppenheimer one. We considered first the XRD structure of hydrated Na vermiculite as reported in ref 16. The space group of vermiculite is monoclinic C2, with cell parameters a = 5.358 Å, b = 9.232 Å, c = 14.97 Å, and β = 96.82°. The simulation box was monoclinic, with lattice parameters a′ = 2 a = 10.716 Å, b = 9.232 Å, c = 14.97 Å, and β = 96.82°, corresponding to two crystallographic cells repeated along a, and included 84 sheet atoms and 4 Na+ cations. The water molecules were 16 and the atoms in total were 136. In other words, the simulated vermiculite had the structural formula Na14H2O[Si3Al1][Mg3]O10(OH)2, corresponding to an asymmetric unit of the crystal, and the simulation box contained four asymmetric units. The structure of the aluminosilicate sheet is given unambiguously in refs 15 and 16, except for the hydrogen atoms forming hydroxyls linked to the octahedral units. The positions of these atoms along z could be determined by neutron diffraction. The density profile of hydrogens for Na vermiculite, reported in ref 19, is blurred in the z range where they should be located, because of the superposition with hydrogens belonging to water. Their position is reported in ref 20, about the structure of Ca vermiculite, where water is tightly coordinated to the cation and the hydrogens of the hydroxyls become visible. In the simulations their positions were optimized by performing a short CPMD run at room temperature starting from reasonable values, by considering the usual OH distance for the acidic protons in aluminosilicates. Starting from the ordered structure of the interlayer hydrated Na+ proposed by Beyer and Graf von Reichenbach,16 as expected on the basis of the force equilibrium considerations reported in the Introduction, a sudden rearrangement of the interlayer content led to a completely disordered structure, even involving the aluminosilicate sheets. Then, we managed to build a more symmetric starting structure, by supposing that the different Na+ cation positions with partial occupancies correspond, as often occurs in crystallography, to a mosaic of different homogeneous domains separated by disordered interfaces. For the cations the m1 site only was considered. As initial positions of water oxygens we considered the sites 1, 2, and 3 reported in ref 16, which are compatible with the positions of the cations, and for the hydrogens we took into account the usual water molecule geometry. As the new structure included only one kind of cation positions, equal intercation distance was ensured by the translational crystal symmetry. The water molecules resulted again to be distributed

atoms each. The chains would be separated by void spaces. In addition, a special choice of the cation occupancy was made (see also Figures 2 and 3). As the asymmetric unit of the crystal contains only one Na+ cation, the nonprimitive C2 unit cell must contain two cations, distributed in four sites, m1, m2, m1′, m2′, the primed ones generated by a (1/2,1/2,0) translation. In ref 16, it was supposed that m1 and m2′ would be fully occupied, whereas m1′ and m2 would be void. Although the proposed structure is compatible with the diffraction data and the geometrical constraints, there is a problem about the equilibrium of the electrostatic forces. Indeed, the distance between Na+ cations, would not be equal, as the m1−m2′ site distance is about 3.1 Å, while the distance between the cation sites with some of their periodic image is about 5.4 Å (see Figure 3). Electrostatic repulsive forces between cations, which cannot be equilibrated by the much weaker water−water and water−sheet hydrogen bonds (HBs) (see Figures 2 and 3), would immediately drive them to a more symmetrical structure where all the interionic distances are equal and Na(H2O)6 octahedra would be rearranged in consequence. Moreover, hydration would reduce or spread the cation charge, but it is well-known that the charge transfer from hydrated Na+ cations to water is negligible.17,18 In a series of experimental papers,19−21 neutron diffraction studies of the interlayer water in vermiculite coordinated to different cations, including Na+, were reported. The results were given as density profiles ρ(z) along z (or c), the axis perpendicular to the sheets, and, by using H/D substitution, it was possible to separate the hydrogen distribution (of octahedral hydroxyls of the sheets and of water) from the oxygen plus clay distribution, including cations. An interesting feature of ρ(z) is its non-negligible value around the plane at z/ c = 0.5, equidistant from the sheets, both for hydrogens and oxygens, indicating the presence of some water in this plane, which was not detected by X-rays. In addition, different from the data reported in ref 15, no water is visible at z/c = 0.3158 (z = 4.724 Å). In this contribution, vermiculite is studied making use of the first-principles Car−Parrinello Molecular Dynamics (CPMD) simulations.22 Ab initio methods have been applied successfully to the study of hydrated clays, and of several physicochemical properties of these materials and in particular of the nanoconfined water.23−34 A review of the investigations performed until 2005 can be found in ref 25 and some details of previous ab initio studies of clays are reported in the Supporting Information. Given the aforementioned disagreement among experimental measurements of its structure, here we focus our simulations on the Na vermiculite clay containing 4 H2O per Na+. In particular structural and vibrational properties are calculated to obtain a consistent picture of the system, and to validate experimental outcomes.15,16,19,21,35−41

2. MODEL AND CALCULATIONS Ab initio MD simulations were performed by using the Car− Parrinello (CP)22 scheme for propagating the wave functions and the ionic configurations as implemented in the CPMD package.42 The PBE density functional43 was used for the electronic structure calculations, together with Vanderbilt ultrasoft pseudopotentials,44 and a wave function cutoff of 25 Ry (32,79 kJ/mol). The Na, Mg, Al, and Si pseudopotentials were constructed from an ionic configuration with a charge of +0.5 e, and the 2s and 2p semicore electrons of Na and Mg 15584

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

among (H2O)6 octahedra, coordinated to the cations and sharing two oxygen atoms each, but with a distribution more uniform than in the structure proposed in ref 16. The new system was simulated again for 75 ps, after equilibration, and the average temperature was 290 K. The resulting simulated structure at room temperature was still disordered, but, as will be shown, average coordinates and density distributions were compatible with experimental data. A bulk system consisting of 64 water molecules was also simulated for reference, using the same pseudopotentials and wave function cutoff as for vermiculite. The cubic box lattice parameter was set to 12.42962 Å in order to reproduce a density of 0.997 g cm−3. Production run in the microcanonical ensemble spanned 15 ps and the average temperature was 302 K. It should be remarked that no symmetry was imposed for the motion of the atoms, but a special procedure was used to verify if, on average, the experimental crystal symmetry was maintained (see the Supporting Information). Average OH distances and HOH angles of water molecules were computed directly via their distribution functions. Vibrational spectra were derived from Fourier transforms of the velocity autocorrelation function of the atomic nuclei. This method permits the evaluation of the vibrational properties of selected groups of atoms separately, such as the aluminosilicate sheets, the cations, or the water molecules. More details and references are given in the Supporting Information.

Table 1. Average Crystallographic Coordinates of Vermiculite: Atoms of the Clay Sheet x/a

atom Si, Al

Si

O1

O2

O3

O4

O5

O6

H

3. RESULTS AND DISCUSSION a. Structure of Vermiculite Sheet and of Adsorbed Water. The computed average coordinates of the atoms belonging to the phyllosilicate sheet are well reproduced (see Table 1). Experimental values in refs 15 and 16 are given without errors as usual in crystallographic papers. This issue is briefly discussed in the Supporting Information. The initial configuration of the sheet was taken from ref 15 and the coordinates of all atoms remained close to the initial values, with oscillations within 0.8 Å and averages differing at most by 0.1 Å from the experimental values. All the coordinate distributions are narrow, Gaussian-like, and unimodal, indicating that the symmetry of the simulated structure is the same as the experimental one (see the Supporting Information). The zcoordinates are in agreement also with the neutron-diffraction density profile reported in ref 19 (see Figure 1a). Also the zcoordinate of the octahedral hydroxyls is in good agreement with experimental data.20 Hydroxyl bonds are approximately parallel to z and point toward the silicate tetrahedra (or away from the octahedral cations). The structure of the interlayer (or intersheet) hydrated Na+ cations is disordered, showing atomic positions with mean square displacements (MSDs) of about 1, 2, and 2.6 Å2 for the cations, water oxygens, and water hydrogens, respectively, for a simulation lasting 75 ps. The plots of the MSDs vs time are not bounded (they show a linear trend) and could suggest the presence of some diffusion (with diffusion coefficient of water of the order of 10−12 m2 s−1), but the simulation is exceedingly too short to ascertain any diffusive behavior and to yield a reliable value. It is also possible that the motion of water and cations is localized around stable positions. We shall consider performing long classical MD simulations in the future to verify these issues, using a reliable force field able to reproduce the CPMD results.

MgI

MgII

MgIII

O3−H (Å) (avg)

exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptlc calcd exptla exptlb calcd exptla exptlb calcd exptla exptlb calcd exptlc calcd

y/b

z/c

0.3974 0.3949 0.3949 0.3858 0.3949 0.3964 0.3590 0.3560 0.3553 0.3750 0.3560 0.3613 0.3580 0.3588 0.3556 0.1421 0.1475 0.1584 0.1471 0.1474 0.1477 0.4407 0.4275 0.4232

0.9980 0.9992 0.9982 0.3303 0.3345 0.3322 0.0016 −0.0030 −0.0021 0.3372 0.3326 0.3320 0.6713 0.6648 0.6671 0.4030 0.4063 0.4069 0.9297 0.9275 0.9205 0.1649 0.1669 0.1697

0.3859 0 0 0.0013 0 0 −0.0018 0 0 −0.0018 0.95 0.93

0.6603 0.9971 0.9971 0.9992 0.3330 0.3324 0.3338 0.6665 0.6648 0.6640

0.1838 0.1830 0.1873 0.1838 0.1830 0.1890 0.0753 0.0717 0.0734 0.0715 0.0717 0.0788 0.0721 0.0689 0.0715 0.2259 0.2221 0.2299 0.2244 0.2221 0.2303 0.2257 0.2212 0.2303 0.137 0.1348 0 0 0.0016 0 0 0.0015 0 0 −0.0014

a

Experimental values are taken from ref 15. bExperimental values are taken from ref 16. cExperimental values are taken from ref 20.

In Figure 1a the experimental and computed density profiles, ρ(z)s, are plotted along with their integrals giving the cumulative densities [the integral of the ρ(z)s]. To compare the computed density profile with the neutron diffraction ρ(z), the contribution of the different chemical species was weighted according to their scattering intensity derived from the data reported in Table 3 of ref 19. Whereas the integral of the ρ(z) of the sheet atoms, water oxygens (Ow), and Na+ cations reaches the same value as the experimental one, the integral of the experimental hydrogen atoms ρ(z) is smaller than the computed value because, as discussed in ref 19, possibly due to an incomplete H/D exchange, the resulting number of hydrogen atoms per structural formula it is only 7.86, whereas in our model it is 10, including 4 water molecules and 2 hydroxyls. We remark that the structural formula reported in ref 16 contains 1.33 hydroxyls and small quantities of Ti, Al, and Fe in the octahedra bearing positive charges, but to reach electrical neutrality without these extra elements 2 hydroxyls are 15585

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

required. The computed density profile of the interlayer content is slightly asymmetric, and thus, in order to make the comparison with experimental data easier, it was symmetrized according to the formula: ρsymmetrized (z) =

1 [ρ(z) + ρ(c − z)] 2

(1)

where c = 14.97 Å is the cell side along z. The agreement with the neutron-diffraction ρ(z) reported in ref 19 is reasonably good including, in the interlayer region, a non-negligible density both of water hydrogens and oxygens around the plane at z/c = 0.5, or z = 7.48 Å (Na+ and Ow contribute to the Na+ + Ow ρ(z) with weights of 0.34 and 0.58, respectively, see Table 3 of ref 19). The experimental and computed cumulative densities of the sheet atoms are in good agreement for z < 3.2 Å, but the experimental one for larger z becomes larger, due to a larger oxygen density peak around z = 3.5 Å. In ref 19 this oxygen excess is ascribed to water adsorbed in the hexagonal cavities present in the basal plane of the sheets, evidenced also by the peaks of the hydrogen ρ(z) around 2.75 and 3.54 Å. In our simulations, as well in other ab initio studies,23−34 water did not reach such positions, but one cannot exclude that it would occur in much longer simulations. In the X-ray structures water was not detected either in these positions or around the plane at z/c = 0.5. The density peak of water oxygens at 6.12 Å (z/c = 0.409) and at 6.06 Å (z/c = 0.405) for the calculated and experimental ρ(z), respectively, closely reproduces the z-coordinate of water oxygens reported in refs 15 and 16 (z/c = 0.4026−0.4045). Another water hydrogen sharp peak is visible in the interlayer region at 5.13 ad at 5.10 Å for the calculated and experimental ρ(z), respectively. Their distance from the basal oxygens of the sheet is 1.7−1.8 Å, a value typical of HBs, whose presence will be discussed below. Of the four Na+ cations, in our simulations two remain near the plane at 7.48 Å (z/c = 0.5), which is the z-coordinate found in the X-ray studies, but the other two are located near the basal oxygens of the sheet (see Figures 1b and 2), at the expected distance of about 2.4 Å, corresponding to the typical value for Na+−Ow, Ow belonging to the first hydration shell. This finding is compatible on the one hand with the total occupancy of Na+ sites reported in the X-ray diffraction experimental papers (0.24 and 0.476 in refs 15 and 16, respectively) and, on the other hand, with the considerations reported in ref 19: “From our experiments we cannot deduce that the interlayer cations lie midway between the clay layers [...]. Rather, the presence of positively charged hydrogen atoms midway between the sheets suggests that some interlayer Na+ lie closer to the clay surface.” X-ray structures report an average value z/c = 0.5, which is still compatible with our estimate z/c = 0.50 ± 0.15 (see Table 2). This average results from the multimodal distribution shown in Figure 1b, which is again compatible with X-ray data, yielding a total occupancy of Na+ cations less than 0.5 at z/c = 0.5. The computed coordinates of the Na+ cations in the xy plane do not obey the crystal symmetry (see also Figures 2 and 3), but for three of the four cations the average positions broadly correspond to the m1 site, and for the other to the m2 site (see Table 2). The problem of having two different values for the interionic distance, as it appears from the Na+−Na+ radial distribution function (RDF) reported in Figure 1c, is explained qualitatively in Figure 3. When at least one water molecule is interposed between two cations, the distance is around 5.5 Å. When no water is interposed the interionic distance is around 3.7 Å. These distances correspond approximately to the

Figure 1. System density profiles and distributions of Na+ cations: (a) ρ(z) (left scale). Black solid line: computed for sheet plus water oxygens and Na cations; red dotted line: the same from neutron diffraction experiments;19 green dashed-dotted line: computed for hydrogens; blue dashed line: the same from neutron diffraction experiments.19 The same symbols hold for the cumulative densities [the integral of the ρ(z)’s] (right scale). (b) Blue solid line (left scale): computed ρNa(z); red dotted line (right scale): cumulative density of Na+ cations. (c) Blue solid line (left scale): Na+−Na+ radial distribution function; red dotted line (right scale): neighbor cation number n(r).

distance between m1 or m2 sites and their periodic images (about 5.4 Å), and to the m1−m2′ site distance (about 3.1 Å), respectively, found in X-ray diffraction experiments.15,16 Notice that in the coordinate optimization procedure starting from diffraction data the space symmetry of the crystalline sheets is assumed, and that the electron density map of the interlayer content is symmetrized accordingly. It is possible that much longer simulations would lead to a more symmetric density map, but the presence of two different Na+−Na+ distance values emerges already from our relatively short trajectory. The current neighbor cation number n(r) plotted in Figure 1c (right scale) shows that the shorter distance is about twice as frequent as the larger one, because the minimum of the first peak corresponds to n(r) ≈ 2. As for the positions of water molecules in the xy plane, the computed average values are still compatible with the experimental ones within their error, given by the root-meansquare displacement over the MD trajectory (see Table 2). Experimental values are also characterized by an overall partial 15586

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

for the different cations and, as verified, in time (see the Supporting Information, Figure S2). Therefore, the RDF plotted in Figure 4 (top panel) carries the correct information about the distribution of Na+−Ow distances, but as the system is strongly anisotropic the current cumulative coordination number is not meaningful, as discussed in ref 52. Fortunately, the cations are only four, and a more detailed analysis can be performed. In the bottom panel of Figure 4 the distributions of coordination number for each Na+ cation is plotted. The coordination number n(Na) ranges from 1 to 7 (see Figure S3 in the Supporting Information), but the distributions are bimodal, except for Na1, which yet is asymmetric. Maxima at 2 and 3 molecules per Na+ are characteristic of partially hydrated cations, when they are close to the sheets. The maxima at 5 molecules per Na+ are visible in all the distributions, and correspond to hydrated cations in the interlayer. Only the Na3 cation shows a significant number of configurations with n(Na) = 6 and, along with the Na2 cation, a few configurations with n(Na) = 7 (see the Supporting Information, Figure S2). Therefore, the presence of stable Na (H2O)6 octahedra proposed in ref 16 is not supported by our calculations. These findings can be compared with those of recent accurate ab initio calculations of the coordination number of Na+. Examples are salt water in carbon nanotubes,53 where n(Na) ranges from 3 to 8, with an average of 4.9 (in reference bulk 5.2). In bulk water solution54,55 the average of n(Na) varies from 4.9 to 6.1, depending on the DFT functional, but with monomodal distributions ranging from 4 to 8.54 Therefore, it seems that small coordination numbers can be related to confinement, independently of the confining system. In Figure 5 the RDF for the distances between water hydrogens and the oxygen atoms of the aluminosilicate sheets (OS) are plotted. The first peaks of OS−Ow (Figure S4 in the Supporting Information) and OS−Hw RDFs fall at 2.72 and 1.72 Å, respectively, indicating the presence of a HB between water molecules and the oxygen atoms of the sheets. Both values are slightly larger than the ones for Ow−Ow and Ow−Hw RDFs (2.68 and 1.66 Å, respectively, see Figure 6), probably indicating a HB slightly weaker than the water−water one. Finally, the Ow−Ow RDF (Figure 6) shows a first neighbors distance larger than the one evaluated with the same model for bulk water (red and shaded in the figures), reflecting the lower density. Accordingly, in the Ow−Hw RDF the first peak is shifted to larger distances. Due to the interactions with the cations and the hydrogen bonds with the oxygens of the sheets a number of water−water HBs are lacking, and in these situations the Ow−Ow distance becomes larger, as it appears from the second anomalous peak in the Ow−Ow RDF at about 3.3 Å. It is interesting to compare our results with those of refs 26 and 27, where hydrated smectite was studied with use of the DFT ab initio method (details are reported in the Supporting Information). The considered system was similar to that considered in this work, the main difference consisting of a layer charge of −0.5 e per unit cell, thus needing only one Na+ cation per simulation box, corresponding to 2 crystallographic cells and containing 16 water molecules. These molecules were sufficient to fill the interlayer space, because they showed weaker HBs of water with the sheets and a liquid-like HB network except for the few coordinated to the cation. The Na+ cation preferred position was close to the sheet, but also the midpoint of the interlayer was stabilized by an energy barrier of ∼1.2kT between the two minima.

Figure 2. Structure of vermiculite as seen from the a axis. Top: A snapshot of the simulated structure (the simulation box is periodically repeated along a and b) showing that some cations (blue) can approach the sheet surface (golden spheres). Bottom: The structure according to ref 16 (same symbols as top panel).

occupancy (at most 0.65), an indication of a remarkable disorder. In Figure 3 the snapshots of the interlayer content in the xy plane are compared with the structure of the interlayer content according to X-ray diffraction experiments,16 above and below the plane at z/c = 0.5. In both simulated and experimental structures water molecules are arranged in distinguishable stripes along the a axis. However, in the simulated structure the water molecule stripes above and below the half intersheet plane seem to be shifted along the b axis and their arrangement is not regular, whereas in the experimental structure the stripes above and below the same plane are more regular and nearly superimposed (see also Figure 2). Clearly, in all cases there is room to accommodate more molecules. In Table 2 we report also the computed average O−H distance and HOH angle, which are 0.997 ± 0.006 Å and 104.5 ± 1.0°, respectively. These values are similar to those found for water coordinated to Na+ in accurate structures (obtained mainly by neutron scattering experiments) of hydrated aluminosilicates, and of zeolites in particular, at least for OH distance, whereas for the HOH angle one would expect a somewhat larger value.49−51 Na+−Ow RDF is reported in Figure 4 (top panel). In spite of the disordered arrangement of water, the molecules are coordinated to the cations as usual, with an average Na+−Ow distance of 2.29 Å, close to the average of the values reported in ref 15, 2.33 Å (see also Table 2). Na+−Hw RDF (see the Supporting Information, Figure S1) has a maximum at 2.89 Å, indicating that hydrogen atoms point mainly out of the cation, as expected. In Figures 2 and 3 it appears that the number of water molecules coordinated to Na+ cations is largely variable 15587

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

Figure 3. The interlayer content as seen from the c axis. Top left: A snapshot of simulated structure. Top right: Structure according to ref 16, below the plane at z/c = 0.5. Bottom: The same respectively above the plane at z/c = 0.5. The blue dots represent the Na+ cations.

dehydrated samples were used,37,38 but they can still be compared with our simulated spectra of the aluminosilicate sheet, which is only slightly affected by hydration, because the interactions of water with the sheet are weak, as seen from the comparison between the experimental data reported in Table 3. In experimental spectra the hydroxyl sharp stretching bands can be distinguished from the underlying broad stretching bands of water, as evidenced by the spectra taken under progressive dehydration.37 To derive the vibrational spectra of the aluminosilicate sheet, the Fourier analysis was performed by considering only the velocities of the atoms belonging to the sheet itself (including Na+ cations), and, for water, only water molecule atoms were included in the calculation. In Table 3 experimental and computed frequencies of the aluminosilicate sheet are compared. The agreement seems to be satisfactory. The frequencies at 114 and 167 cm−1 are to be ascribed to the vibrations of Na+ cations, as verified by performing the Fourier transform of the cation velocities only. Frequencies in the range 180−1200 cm −1 correspond to the vibrations of the aluminosilicate sheet skeleton, while those in the range 1400−3800 cm−1 correspond to the bending (1400−1700 cm−1) and to the stretching (2900−3800 cm−1) modes of octahedral OH bonds. A detailed analysis of the vibrational mode assignment is reported in ref 38. Overall, the standard error of the computed frequencies, with respect to the average experimental values, is of the order of 4−5% for all the

Table 2. Crystallographic Coordinates, Distances and Angles: Interlayer Content, Na Ion, and Water Oxygensa x/a

atom Na1 Na2 Ow1 Ow2 Ow3 Na−O (Å) (avg) Ow−H (Å) H−Ow−H (deg)

exptlb calcd exptlb calcd exptlb calcd exptlb calcd exptlb calcd exptlb

0.5000 0.66 ± 0.5000 0.39 ± 0.1546 0.17 ± 0.1561 0.19 ± 0.1546 0.17 ± 2.33

calcd calcd calcd

2.29 0.997 ± 0.007 104.4 ± 1.2

0.28 0.10 0.21 0.11 0.07

y/b 1.0000 0.91 ± 0.3370 0.39 ± 0.0169 0.09 ± 0.3121 0.34 ± 0.6667 0.69 ±

z/c 0.20 0.05 0.08 0.03 0.11

0.5000 0.50 ± 0.15 0.5000 0.50 ± 0.15 0.4045 0.43 ± 0.04 0.4045 0.43 ± 0.04 0.4045 0.407 ± 0.003

a

Errors of computed coordinates correspond to the square root of MSDs. A comment about the errors of experimental coordinates can be found in the Supporting Information. bExperimental values are taken from ref 15.

b. Simulated Vibrational Spectra. The experimental vibrational spectra of vermiculite were recorded with use of IR35,38−41 and Raman36−38 spectroscopies. In some cases 15588

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

Figure 4. Top: Computed Na+−Ow radial distribution function (RDF). The red arrow shows the average of the experimental Na+−Ow distances. Bottom: Distributions of coordination number for each Na+ cation. Symbols represent histogram values, with unitary bins. Lines are to help guide the eye.

Figure 6. Computed water−water RDFs in vermiculite. Top: Oxygen−hydrogen RDF. Bottom: Oxygen−oxygen RDF. The lines including the shaded areas are the RDFs evaluated for bulk water.

Table 3. Experimental and Computed Frequencies (cm−1) of Aluminosilicate Sheeta exptl IR (ref 35)

exptl Raman (ref 36)

exptl Raman (ref 37)

105 164

exptl IR (ref 38)

Figure 5. Computed OS (silicate sheet oxygens)−Hw RDF. 676 723 812 957 995 1095

3220 3550

exptl IR (ref 40)

128 151

367 429 458 508

frequencies up to 1650 cm−1, but it lowers to about 2% for the stretching vibrations of octahedral OH bonds. In Figure 7 we report the vibrational spectra of adsorbed water in low-frequency, bending and stretching regions, respectively. The computed spectra are compared with the corresponding ones evaluated for bulk water at room temperature. Overall, the spectrum of adsorbed water is similar to that of bulk water (reference CPMD simulations). However, as expected from the general trends for confined water coordinated to metallic cations (for instance in zeolites), low-frequency spectra of adsorbed water, corresponding to librational modes are overall slightly red-shifted with respect to those of bulk water, because in the adsorbed water some HBs are weaker or lacking. In addition, one observes a slight blue-shift of the bending frequencies (from about 1550 to 1590 cm−1), due to the interaction with the cations, and the presence of two main peaks in the stretching band, the more blue-shifted being

exptl Raman (ref 39)

3640 3769

187 278 310 354 434 532 552 634 678 741 960 1033 1083

456

686 758 978 998

1643

1422 1622

3207 3390 3561 3711

3427 3558 3670

calcd (this work) (114) 167 186 256 (300) 378 418 465 (506) (550) 630 664 752 880 935 971 1470 1590 3000 3238 3380 3591 3627

a

For computed frequencies, values in parentheses are shoulders or weak peaks.

15589

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

sufficiently long simulation time (about 75 ps) resulted in a good reproduction both of the structure and of vibrational properties of vermiculite aluminosilicate sheet. We focused on the analysis of the structure of the interlayer content (Na cations and water molecules) because experimental results, obtained by X-ray and neutron diffraction at first sight, do not agree. The former technique is biased by the crystalline symmetry of the aluminosilicate sheet, and the optimized coordinates are affected by partial occupancies, indicating that some positions or distributions are missing. The density profiles perpendicular to the sheets (for hydrogens and for the other elements) derived from neutron diffraction yields complementary information, such as the presence of both water hydrogens, and Na+ cations and/or water oxygens around the plane at z/c = 0.5 parallel to the sheets and the presence of HBs between the adsorbed water and the basal oxygens of the aluminosilicate sheets. From the simulation results the interlayer content appears to be remarkably disordered, but the average coordinates are compatible with the X-ray diffraction data, taking into account also the partial occupancies of the sites, and the computed density profile is similar to that reported in the neutron diffraction studies. Although the water molecules are coordinated to the cations in the usual way, they form incomplete solvation shells around the cations and account for the presence of different values of Na+−Na+ distances, as found in the X-ray studies, corresponding to the presence, or not, of water molecules interposed between the cations. In addition, the lack of some HBs between water molecules and the presence of HBs with the basal oxygens of the aluminosilicate sheets result in weakened interaction between some water molecule pairs, which is evidenced by an anomalous peak in the O−O RDF at about 3.3 Å. In addition, it appears that there is room to accommodate more water molecules in the water layers. We have shown, therefore, that, besides its partial disagreement with neutron diffraction data, the ordered structure proposed in ref 16 hardly would correspond to the real one. Indeed, not only do the electrostatic forces between cations not balance, but also the presence of stable Na(H2O)6 octahedra is questionable because of the scattered positions of the Na+ cations, some of them being partially coordinated to the basal oxygens of the sheet. The comparison with the ab initio calculations of hydrated smectite clay23−34 can cast some light on the behavior of water in clays in general. In smectite the sheet charge is smaller and for a system of the same size as ours only one Na+ cation is present along with 16 water molecules, and the results of the calculations show that it is just the sheet charge that causes the different interlayer content structure. The larger charge in vermiculite entails stronger HBs between water and sheet surface, and a 4-fold number of Na+ cations, which in turn coordinate more water molecules, thus resulting in weaker water−water HBs, smaller water mobility, but also more available space and disorder. On the contrary, in smectite the adsorbed water is liquid-like, as there are less constraints stemming from the HBs with the sheet and from the cation coordination, also because the preferred position of the Na+ cation is near the sheet. There is an open question about the water adsorbed in the hexagonal cavities present in the basal plane of the sheets, resulting from the analysis of the neutron diffraction data.19 As remarked above, neither in this study nor in previous ab initio

Figure 7. Simulated IR spectra of water in vermiculite (blue dotted line) and reference CPMD simulations for 64 water molecules (black solid line). Top: Low-frequency region, corresponding to translations and librations of the water molecules. Middle: The water bending mode region. Bottom: The water stretching mode region.

caused by the lack of some HBs between water molecules. In Table 4 the computed and experimental frequencies for Table 4. Experimental and Computed Frequencies (cm−1) of Adsorbed Water exptl Raman (ref 37)

3450

exptl IR (ref 40)

exptl IR (ref 41)

1640

1650

3400

3210 3450

calcd (this work) 1590 3019 3198 3372 3461

bending and stretching bands are collected. The relative error with respect to the average of the experimental values is 4.4% and 2.0% for bending and stretching bands, respectively, similar to that of the aluminosilicate sheet. Therefore, the stretching vibration frequencies of OH bonds is reproduced with a noticeable accuracy.

4. CONCLUSIONS The above-reported CPMD simulations of vermiculite with the experimental water content (eight H2O per u.c., or four per Na+ cation) show that the choice of using the PBE density functional43 for the electronic structure calculations, together with Vanderbilt ultrasoft pseudopotentials,44 as well as of a 15590

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

(15) Slade, P. G.; Stone, P. A.; Radoslovich, E. W. Clays Clay Miner. 1985, 33, 51−61. (16) Beyer, J.; Graf von Reichenbach, H. Clay Miner. 2002, 37, 157− 168. (17) Collins, K. D.; Neilson, G. W.; Enderby, J. E. Biophys. Chem. 2007, 128, 95−104. (18) Collins, K. D. Biophys. J. 1997, 72, 65−76. (19) Skipper, N. T.; Soper, A. K.; McConnell, J. D.C. J. Chem. Phys. 1991, 94, 5751−9574. (20) Skipper, N. T.; Soper, A. K.; Smalley, M. V. J. Chem. Phys. 1994, 98, 942−945. (21) Sposito, G.; Skipper, N. T.; Sutton, R.; Park, S.; Soper, A. K.; Greathouse, J. A. J. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 3358−3364. (22) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471−2474. (23) Boulet, B.; Greenwell, H. C.; Stackhouse, S.; Coveney, P. V. J. Mol. Struct.: THEOCHEM 2006, 762, 33−48. (24) Stackhouse, S.; Coveney, P. V.; Sandré, E. J. Am. Chem. Soc. 2001, 123, 11764−11764. (25) Tunega, D.; Benco, L.; Haberhauser, G.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. B 2002, 106, 11515−11525. (26) Boek, E. S.; Sprik, M. J. Phys. Chem. B 2003, 107, 3251−3256. (27) Refson, K.; Park, S.-H.; Sposito, G. J. Phys. Chem. B 2003, 107, 13376−13383. (28) Stackhouse, S.; Coveney, P. V.; Benoit, D. M. J. Phys. Chem. B 2004, 108, 9685−9694. (29) Chatterjee, A.; Ebina, T.; Onodera, Y.; Mizukami, F. J. Chem. Phys. 2004, 120, 3414−3424. (30) Suter, J. L.; Boek, E. S.; Sprik, M. J. Phys. Chem. C 2008, 112, 18832−18839. (31) Clausen, P.; Andreoni, W.; Curioni, A.; Hughes, E.; Plummer, C. J. G. J. Phys. Chem. C 2009, 113, 12293−12300. (32) Clausen, P.; Andreoni, W.; Curioni, A.; Hughes, E.; Plummer, C. J. G. J. Phys. Chem. C 2009, 113, 15218−15225. (33) Liu, X.; Meijer, E. J.; Lu, X.; Wang, R. Clays Clay Miner. 2010, 58, 89−96. (34) Churakov, S. V.; Kosakowski, G. Philos. Mag. 2010, 90, 2459− 2474. (35) Serratosa, J. M.; Bradley, W. F. J. Phys. Chem. 1958, 62, 1164− 1167. (36) Wada, N.; Suzuki, M.; Hines, D. R.; Koga, K.; Nishihara, H. J. Mater. Res. 1987, 2, 864−870. (37) Wada, N.; Kamitakahara, W. A. Phys. Rev. B 1991, 43, 2391− 2397. (38) Arab, M.; Bougeard, D.; Smirnov, K. S. Phys. Chem. Chem. Phys. 2002, 4, 1957−1963. (39) Xie, Y.; Wang, A. J. Compos. Mater. 2009, 43, 2401−2417. (40) Tomanec, R.; Popov, S.; Vučinič, D.; Lazič, P. Fiz. Probl. Mieralurgii 1997, 31, 247−254. (41) Muiambo, H. F.; Focke, W. W.; Atanasova, M.; van derWesthuizen, I.; Tiedt, L. R. Appl. Clay Sci. 2010, 50, 51−57. (42) CPMD, Version 3.12; copyright IBM Corp., 1990−2009, copyright MPI für Festkörperforschung Stuttgart, 1997−2001; www. cpmd.org. (43) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (44) Vanderbilt, D. Phys. Rev. B 1990, 41, 7892−7895. (45) Van de Vondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Phys. Chem. B 2004, 108, 12990−12998. (46) Nosé, S. Mol. Phys. 1984, 52, 255−268. (47) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (48) Martyna, G. J.; Tuckerman, M. E.; Klein, M. L. J. Chem. Phys. 1992, 97, 2635−2643. (49) Artioli, G.; Smith, J. V.; Kvick, Å. Acta Crystallogr., Sect. C 1984, 40, 1658−1662. (50) Peacor, D. R. Am. Mineral. 1973, 58, 676−680. (51) Torrie, B. H.; Brown, I. D.; Petch, H. E. Can. J. Phys. 1981, 74, 229−240. (52) Demontis, P.; Gulín González, J.; Suffritti, G. B. J. Phys. Chem. C 2012, 116, 11100−11109.

works was water found in these positions, but possibly calculations did not consider them and molecular dynamics simulations were too short, at most some tenths of picoseconds, so that much longer simulations, or simulations run at high temperature, could drive the water molecule to reach them. On the experimental side, it is possible that these kind of positions require overcoming an energy barrier, so that their population could depend on the thermal history of the samples. Work is in progress to extend the simulations to vermiculite containing more water molecules. Preliminary results show that more ordered arrangement of water with interesting features is obtained.



ASSOCIATED CONTENT

S Supporting Information *

More details and references. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*Fax: +39-079-229559. Tel: +39-079-229552. E-mail: pino@ uniss.it. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This research is supported by the Italian Ministero dell’Istruzione, dell’Università, e della Ricerca (MIUR), PRIN funding, by Università degli studi di Sassari and by Istituto Nazionale per la Scienza e Tecnologia dei Materiali (INSTM), which are acknowledged. The “Consorzio COSMOLAB” is also acknowledged for the resources provided within the Cybersar Project, as well as the CASPUR project std12-031 “Ab initio simulations of water confined in vermiculite clay”. The authors are thankful to Prof. Bernd Meyer for having kindly provided Vanderbilt ultrasoft pseudopotentials for this study and for useful discussion.



REFERENCES

(1) Sposito, G. The Chemistry of Soils; Oxford University Press: New York, NY, 1989. (2) Hydrous Phyllosilcates; Bailey, S. W., Ed.; Mineralogical Society of America: Washington, D.C., 1988. (3) Moore, D. M.; Reynolds, R. C. X-Ray Diffraction and Identification and Analysis of Clay Minerals; Oxford University Press: New York, NY, 1997. (4) Borchardt, G. Minerals in Soil Environments; Dixon, J. B., Weed, S. B., Eds; Soil Science Society of America: Madison, WI, 1989; p 675. (5) Sposito, G. The Surface Chemistry of Soils; Oxford University Press: New York, 1984. (6) Karaborni, S.; Smit, B.; Heidug, W.; Urai, J.; van Oort, E. Science 1996, 271, 1102−1104. (7) Mishima, O.; Stanley, H. E. Nature 1998, 396, 329−335. (8) Bergman, R.; Swenson, J. Nature 2000, 403, 283−286. (9) Bergman, R.; Swenson, J.; Börjesson, L.; Jacobsson, P. J. Chem. Phys. 2000, 113, 357−363. (10) Swenson, J.; Jansson, H.; Howells, W. S. J. Chem. Phys. 2000, 113, 2873−2879. (11) Swenson, J.; Bergman, R.; Longeville, S. J. Chem. Phys. 2001, 115, 11299−11305. (12) Jansson, H.; Swenson, J. Eur. Phys. J. E: Soft Matter Biol. Phys. 2003, 12, S51−S54. (13) Swenson, J.; Jansson, H.; Howells, W. S.; Longeville, S. J. Chem. Phys. 2005, 122, 084505. (14) Swenson, J.; Teixeira, J. J. Chem. Phys. 2010, 132, 014508. 15591

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592

The Journal of Physical Chemistry C

Article

(53) Kulik, H. J.; Schwegler, E.; Galli, G. J. Phys. Chem. Lett. 2012, 3, 2653−2658. (54) Bankura, A.; Carnevale, V.; Klein, M. L. J. Chem. Phys. 2013, 138, 014501. (55) Rowley, C. N.; Roux, B. J. Chem. Theory Comput. 2012, 8, 3526−3535.

15592

dx.doi.org/10.1021/jp4007944 | J. Phys. Chem. C 2013, 117, 15583−15592