- © 2015 Clay Mineral Society
Atomistic simulations of 2:1 clay minerals based on parameterized forcefields have been applied successfully to provide a detailed description of the interfacial structure and dynamics of basal planes and interlayers, but have made limited progress in exploring the edge surfaces of these ubiquitous layer-type aluminosilicates. In the present study, molecular dynamics simulations and energy-minimization calculations of the edge surfaces using the fully flexible CLAYFF forcefield are reported. Pyrophyllite provides an ideal prototype for the 2:1 clay-mineral edge surface because it possesses no structural charge, thus rendering the basal planes inert, while crystal-growth theory can be applied to identify two major candidates for the structure of the edge surfaces. Models based on these candidate structures reproduced bulk crystal bond distances accurately when compared to X-ray data and ab initio molecular simulations, and the predicted edge surface bond distances were in agreement with those determined via ab initio simulation. The calculated surface free energy and surface stress led to an accurate prediction of pyrophyllite nanoparticle morphology, while surface excess energies calculated for the edge surfaces were always negative. These results are consistent with the observed pyrophyllite nanoparticle morphology, with the concept of negative interfacial energies, and conditions that may give rise to them including a role in the stabilization of layer-type nanoparticulate minerals. Molecular dynamics simulations of hydrated nanoparticle edge surfaces indicated five reactive surface oxygen sites on the dominant candidate edge, in agreement with a recent model of proton titration data for 2:1 clay minerals. These promising results illustrate the potential for classical mechanical atomistic simulations that explore edge surface phenomena at much greater length- and times-scales than are currently possible with computationally expensive ab initio methods.
The edge surfaces of 2:1 clay minerals are the source of their pH-dependent charge, affecting colloidal stability and rheological properties (Rand et al., 1980; Ramos-Tejada et al., 2001; Tombacz and Szekeres, 2004) and playing important roles as adsorbents in nutrient retention (Sposito, 2008), pollutant attenuation (Bergaya et al., 2006), and organic matter stabilization (Stevenson, 1994; Lal, 2004; Wan et al., 2007). The dissolution of clay mineral nanoparticles has been observed to proceed predominantly at edge surfaces (Bosbach et al., 2000; Bickmore et al., 2001), and the edge surface plane represents a boundary that diffusing solutes must cross to move between interlayer nanopores and adjacent micro- or meso-pores (Rotenberg et al., 2007; Hedstrom and Karnland, 2012).
Despite this acknowledged importance of the 2:1 clay mineral edge surfaces, an understanding of their reactivity continues to present a major research challenge because a consensus does not yet exist as to their molecular structure (Bleam et al., 1993; Elzinga and Sparks, 1999; Bickmore et al., 2001; Morton et al., 2001; Tournassat et al., 2004; Churakov, 2006; Dähn et al., 2006; Bourg et al., 2007; Dähn et al., 2011; Churakov and Dähn, 2012). To be sure, edge surface structure and reactivity have been addressed by ab initio quantum mechanical simulations (Bickmore et al., 2003; Churakov, 2006; Liu et al., 2012; Martins et al., 2014), including interfacial dynamics (Churakov, 2007) and acid-base properties (Bickmore et al., 2003; Churakov, 2007; Tazi et al., 2012; Liu et al., 2013; Liu et al., 2014). But these kinds of simulations are computationally expensive and, therefore, limited in the length- and time-scales they can explore. Moreover, they proceed from pristine surfaces cleaved a priori from the bulk mineral structure for what is inherently an anisotropic and disordered material, thus revealing that edge surfaces undergo considerable relaxation (Bickmore et al., 2003; Churakov, 2006; Tazi et al., 2012), leading to a disordered structure exhibiting five-coordinated Al (Liu et al., 2012; Tazi et al., 2012) along with spontaneous proton transfers among surface OH groups (Churakov, 2007).
Molecular simulations based on classical mechanics offer an approach complementary to ab initio simulations by allowing a computationally inexpensive exploration of much larger length- and time-scales than ab initio simulations can currently achieve. Most molecular mechanics simulations to date, however, have been based upon a rigid clay mineral structure along with an ad hoc choice of how to model the edge surface (Rotenberg et al., 2007; Churakov and Dähn, 2012; Hedstrom and Karnland, 2012). Evaluation of alternative edge surface structures using classical mechanical simulations based upon a fully flexible force field which take full advantage of the large length-and time-scales available could, therefore, provide useful insights linking ab initio simulations to molecular spectroscopy or other experimental methods that have been applied to probe the 2:1 clay mineral edge surface.
In the present study, this possibility is addressed by implementing the fully-flexible CLAYFF forcefield (Cygan et al., 2004) as a starting point for molecular dynamics simulations of the 2:1 clay mineral edge surface. The interatomic potentials in CLAYFF are parameterized using structural and spectroscopic data for a wide variety of compounds, while a flexible single point charge (SPC)-based model (Berendsen et al., 1987) is used to describe water molecules and hydroxyls. The 2:1 clay mineral selected was pyrophyllite [Al2[Si4O10](OH)2]. Pyrophyllite provides an ideal prototype for the 2:1 clay mineral edge surface because its basal planes are inert, while crystal growth theory can be applied readily to identify candidate structures for its edge surfaces. Aluminum ions occupy two of three possible sites in the octahedral sheet, with internal OH groups aligned in the cis-configuration (Bergaya et al., 2006). Silicon ions fill the cation sites in the two tetrahedral sheets. Thus, pyrophyllite possesses no structural charge, so only the edge surface is strongly reactive.
Candidate edge surface structures
From the perspective of crystal-growth theory (Hartman, 1973), the structure of a 2:1 dioctahedral clay mineral can be envisioned as a network of bond chains based on a tetrahedron-octahedron-tetrahedron (TOT) unit. In a seminal paper, White and Zelazny (1988) applied this perspective to identify three bond chains that run parallel to axes defined by the ways in which two tetrahedral apices in a TOT unit can be joined to an octahedron across a center of symmetry. Two of the three bond chains, identified as A and C chains by White and Zelazny (1988), are symmetrically equivalent and henceforth are referred to collectively here as the AC chain. The remaining bond chain is designated the B chain following White and Zelazny (1988).
The molecular structure of the AC chain as viewed along the chain direction identified by White and Zelazny (1988) is provided (Figure 1a). The edge of an AC bond chain (right side of Figure 1a) has two singly-coordinated O bonded to tetrahedral cation sites, one transitional O bonded to both octahedral and tetrahedral cation sites, one singly coordinated O bonded to an octahedral cation site, and one doubly coordinated O bonded to an octahedral cation site. This doubly coordinated O is part of an exposed structural OH group that seems to have been omitted inadvertently in previous studies of the 2:1 clay mineral edge surface (White and Zelazny, 1988; Bickmore et al., 2003; Liu et al., 2012). White and Zelazny (1988) introduced the AC edge mineral structure including this OH group (their figure 3), but left it out of the structure in their subsequent figures (their figure 5). The AC edge surface structure depicted here (Figure 1a) corrects this longstanding apparent oversight. The edge of a B bond chain (Figure 1b) has pairs of singly-coordinated O bonded to tetrahedral and octahedral cation sites.
These two types of bond chain define the stable edge surface structures for 2:1 clay minerals, i.e. they are identified by crystal-growth theory as minimizing the number of potentially reactive edge surface sites per unit area, thus minimizing the edge surface free energy (Hartman, 1973; White and Zelazny, 1988). Note that the alignment of edge surface O differs for the AC and B chains. Those at the B edge lie in a single plane perpendicular to the basal plane (Figure 1b) whereas those at the AC edge are distributed across three planes perpendicular to the basal plane, with one O bonded to a tetrahedral site and the doubly coordinated O lying in the innermost plane; the singly coordinated O bonded to an octahedral site and the transitional O lying in the central plane; and the O bonded to the second tetrahedral cation site lying in the outermost plane (Figure 1a). Note that a recent study comparing ab initio and classical mechanical simulations of the pyrophyllite edge surface (Martins et al., 2014) chose structures on the basis of the Miller indices that, due to the non-orthogonal simulation cell angles, only roughly approximate the edges surfaces defined by White and Zelazny (1988). As a result of this choice, the AC edge surface was in fact wholly neglected while other edge surface structures were included that are not consistent with crystal growth theory (i.e. asymmetric attachment of the TOT unit or non-optimal reactive surface group densities).
Edge surface proton distributions
A number of proton-distribution models have been proposed for 2:1 clay mineral edge O sites (White and Zelazny, 1988; Bleam, 1993; Bleam et al., 1993; Bickmore et al., 2003), but all such models are subject to the constraint that the edge surface must be at the point of zero net proton charge (p.z.n.p.c.; Sposito, 2008). Three different postulated charge distributions for the AC edge surface (depicted in figure 2c–e in Bickmore et al., 2003) and a proton distribution for the B edge surface that is a hybridization of models proposed by White and Zelazny (1988) and Bleam et al. (1993) are presented (Figure 1c–f). The algebraic sum of bond strengths plus anionic charge, calculated in accordance with Pauling’s second rule (Sposito, 2008) and labeling each edge site O, can be used to evaluate the stability of a proposed edge surface structure. The two charge distributions for the AC edge proposed by Bleam et al. (1993) localize equal and opposite charges at the singly coordinated and transition O bonded to edge Al (Figure 1c,d). An additional AC edge charge distribution was proposed by White and Zelazny (1988). This distribution postulates a relaxed >AlOH2 bond, such that the O now possesses no unsaturated valence and the transition O becomes saturated due to a compensating increase in bond strength which this site receives from Al. A similar charge distribution for the B edge surface is proposed herein (Figure 1f), combining the bond-relaxation concept of White and Zelazny (1988) with the singly and doubly protonated O bonded to octahedral cation sites proposed by Bleam et al. (1993) for the B edge surface. According to Pauling’s Rules, the first two charge distributions (Figure 1c,d) should be less stable than the third (Figure 1e). Moreover, ab initio simulations of the pyrophyllite AC edge surface (Churakov, 2007; Liu et al., 2014) predict that the transition O is unlikely to be protonated at the p.z.n.p.c.
Edge surface models
Cleavage of the pyrophyllite structure to create edge surfaces bounded by the AC or B periodic bond chains was implemented using two different pyrophyllite polytypes, 1Tc and 2M. This cleavage preserves the stoichiometry of the clay mineral (Bickmore et al., 2003), after which healing the cleaved surfaces with protons, hydroxide ions, or water molecules can be implemented to satisfy exposed valencies while altering the stoichiometry only by an integral number of water molecules (Bleam et al., 1993). Such pristine edge surfaces cut from a bulk crystal structure and healed through adsorption of water molecules must be at the p.z.n.p.c. (Sposito, 2008).
A model of the AC edge surface (30 Å vacuum slab oriented normal to the ac plane) was created using the atomic coordinates of pyrophyllite-1Tc (Lee and Guggenheim, 1981). A 4 × 4 × 2 (a × b × c) expansion of the unit cell cleaved along the (010) plane corresponds to the AC bond chain (Figure 1a; White and Zelazny, 1988). Under-coordinated Si1 and Al were healed with OH− or H2O, whereas unsaturated Si2–O were healed with protons. For simulations of the hydrated surface, the vacuum slab was packed with water molecules to a density of 1.00 g/cm3 .
The B chain edge surface was created using the atomic coordinates of pyrophyllite-2M (Bleam et al., 1993). A 4 × 4 × 2 (a × b × c) expansion of the unit cell cleaved along the (010) plane results in the creation of an edge surface consistent with the B bond chain (Figure 1b). Under-coordinated Si, Al, and Si–O were healed with protons, OH−, or H2O molecules. For simulations under hydrated conditions, a 30 Å vacuum slab oriented normal to the ac plane was packed with water molecules to a density of 1.00 g/cm3 .
The proton distributions proposed for the AC edge surface (Figures 1c,d,e) can be simulated by the CLAYFF forcefield with minimal modification of its atom-type assignments for bulk pyrophyllite. CLAYFF assigns an atom type to each atom in a molecular structure based on the chemical element the atom represents and its coordination environment, and the parameterization is such that no new atom types or revisions of partial charges are necessary to simulate the AC edge surface models. The atom type assignments for two of the structures (Figures 1c,d), however, do not lead to a neutral B edge surface. On the other hand, the atom type assignments for the third AC edge surface charge distribution (Figure 1e) do produce a neutral edge surface for both kinds of bond chain. For this reason, and their greater stability according to Pauling’s Rules, the AC and B edge surface charge distributions with >AlOH2 bonds possessing zero bond strength (Figure 1e,f) were used exclusively in these simulations.
Two kinds of water molecule were required to heal the cleaved AC and B edge surfaces, neutral H2O molecules the atoms of which are assigned in accordance with the extended simple point charge (SPC/E) water model (Berendsen et al., 1987) and water molecules that heal unsaturated surface sites through dissociative chemisorption. The H and O of a dissociated water molecule were assigned their respective CLAYFF OH atom types [ho and oh (Cygan et al., 2004)]. A dissociated water molecule heals the >Si2–O site with a proton and the >Si1 – site with a hydroxide ion (Figure 1e,f). Charge neutrality of the edge surface was then maintained through reclassification of the O at the >Si2–O site from bridging- to hydroxyl-O (ob → oh).
The aqueous phase was described by the SPC/E water model including harmonic bond stretching and bending terms (Berendsen et al., 1987; Lopez-Lemus et al., 2008). The CLAYFF forcefield parameters (Cygan et al., 2004) include both long- and short-range interactions, as well as the two harmonic bond vibrational terms. Metal cation–O interactions are described by a Lennard-Jones 6–12 potential along with a Coulomb term utilizing partial charges derived by quantum mechanical methods. The general form of the energy expression is thus: (1)
The multiplier on the right side contains e, the proton charge, and εo, the permittivity of vacuum (8.85419 × 10−12 F/m2). The first term following represents the Coulomb potential, qi and qj being the quantum mechanical partial charges on atoms i and j, with rij being the separation distance between atoms i and j. The second term is a Lennard-Jones potential, where Do,ij (kcal/mol) is the depth of the potential well and Ro,ij (Å) is the collision distance for an ij pair interaction (Cygan et al., 2004). Lorentz-Berthelot combination rules (Allen and Tildesley, 1989) were used to calculate these parameters for heterogeneous atom pairs. The third and fourth terms represent harmonic bond stretching and bending energies, respectively, where ro and θo define the equilibrium bond length and angle, k1 and k2 are force constants, and θijk is the i-j-k bond angle. Long-range electrostatic interactions and the attractive contribution to the short-range interactions were evaluated using Ewald summation (Karasawa and Goddard, 1989) with the repulsive contribution truncated at a distance of 9.0 Å. All Ewald summations were determined to an accuracy of 0.0001 kcal/mol.
Pyrophyllite edge surface structures were investigated through geometry optimization at 0 K followed by molecular dynamics (MD) simulations in the microcanonical ensemble (NVE) where the number of atoms (N), volume (V), and energy (E) are held constant consistently with temperature (T) equal to 298.15 K (Figure 2). All edge surfaces were geometry-optimized in vacuo prior to hydration. Energy-minimized, 3-D periodic structures were obtained using a smart minimization approach that combined, in sequence, steepest descent, adjusted basis Newton-Raphson, and quasi-Newton methods. All simulation cells were triclinic, with the angles for the AC surface [α = 91.18°, β = 100.46°, γ = 89.64°] taken from Lee and Guggenheim (1981). The angles for the B edge surface [α = 90.19°, β = 110.39°, γ = 89.84°] were optimized to maintain the ab plane approximately parallel to the pyrophyllite basal plane. A 10 ps simulation in the canonical ensemble (NVT), where N, V, and T were held constant, was performed to bring the system to a temperature of 298.15 K. A direct velocity-scaling algorithm was used to maintain temperature within ±10.0 K of the target temperature (Accelrys, 2008). Subsequently, a 100 ps simulation in the isobaric-isothermal ensemble (NPT; P = 1.0 atm and T = 298.15 K); where N, the pressure (P), and T were held constant, was completed to optimize the cell volume. These pre-equilibration simulations were used to determine the volume to be used in the NVE simulations, which were conducted for 50 ps to establish equilibrium before the production simulation of 100 ps. A time step of 0.5 fs was used, with the system configuration recorded every 0.1 ps. After this sequence of simulations to establish a low-energy state was completed, the bulk aqueous phase was removed and a final geometry optimization at 0 K was performed to acquire structural data. Representative frames of the hydrated atomistic models from the NVE-MD are provided (Figure 2). All simulations were performed on the computing cluster at the Molecular Graphics and Computation Facility (MGCF), University of California at Berkeley using the Forcite Plus module in Materials Studio 5.5 (Accelrys Inc., San Diego, California, USA).
Surface energies and radial distribution functions
The surface free energy; surface free energy and stress; and surface excess energy were calculated for both the edge and the basal surfaces. The AC and B edge surface energies were calculated from the expanded pyrophyllite unit cells described above. The basal surface energy calculations used a 4 × 4 × 4 (a × b × c) expansion of the triclinic unit cell to create sufficient bulk to separate the surfaces created. The surface free energy, , defined as the work necessary to form a unit surface area, was calculated according to the expression: (2)where Ebulk is the energy of the geometry-optimized bulk crystal; Epore is the energy required to cleave the geometry-optimized bulk crystal and separate the resulting two interfaces by a given distance (pore width); and A is the area of the surface created. For solids, the surface free energy and surface stress are not equal. The surface stress, Eτ, is the work associated with the reconstruction or stretching of a unit surface area (Gibbs, 1931; Adamson and Gast, 1997). Direct experimental measurement of Eτ is difficult (Adamson and Gast, 1997). Thus, the surface free energy and stress, , are considered together according to the expression: (3)where is the energy of the cleaved crystal with geometry-optimized atomic positions. The calculated values were used in a Wulff construction (Wulff, 1901) to predict nanoparticle morphology using the open-source program Wulffman (Roosen et al., 1998).
The surface excess energy, Eγ, includes the energy contributions of chemisorbed water molecules that heal the broken bonds of a cleaved edge surface. The surface excess energy is calculated according to the expression: (4)where is the energy of the cleaved, relaxed, and healed structure; is the energy of the chemisorption reactions necessary to heal one bond chain unit (i.e. >H2O, >OH, and O>H), with n being the number of edge bond chain units. This term is the sum of the in vacuo energy of the SPC/E water molecule and the energy of the dissociative adsorption of a water molecule to the edge Si. Hess’s law (equation 5) was used to calculate this dissociative adsorption energy in the manner of de Leeuw (2008). The lattice energies of SiO2(s) and Si(OH)4, (Elatt[SiO2] and Elatt[Si(OH)4]) calculated with CLAYFF and the experimental heat of formation for silicic acid (Ramberg, 1954; Greenberg, 1957) were used to complete the energy cycle and quantify the energy of dissociative sorption of a H2O molecule to the edge Si (Ediss [H2O]); the energy of the chemisorption reactions to heal one bond chain unit, , was equal to −435.1 kcal/mol. (5)
The radial distribution function (RDF), gαβ(r) for an edge site-water O, was calculated conventionally as (Enderby and Neilson, 1981; Sposito et al., 1999): (6)where the gαβ(r) gives the variation in atomic density of β as a function of its radial distance, r, from atom species α; V is the volume of the simulation cell; Nβ is the total number of β atom species; and dnαβ is the average number of β atom species in a radial shell with thickness dr at distance r from atom species α. The first minimum of the RDF, ρ, defines the radial boundary of the first coordination shell centered on atomic species α (Enderby and Neilson, 1981; Sposito et al., 1999). The development of an ordered solvent shell surrounding edge O sites and the radial distance to this solvent shell were used as a proxy for edge O reactivity, as it is not unreasonable to equate solvent access with O reactivity. On integration, the RDF gives the coordination number, CNαβ(ρ), the average number of species β found within the limits of integration [0, ρ] for a spherical shell centered at α (Allen and Tildesley, 1989): (7)
Bond lengths in the pyrophyllite structure, geometry-optimized in vacuo, are compared with results from experiment (Lee and Guggenheim, 1981) and simulations (Refson et al., 2003; Cygan et al., 2004; Churakov, 2006) (Table 1a). The bond distances, measured on four periodic bond chains in the center of a layer, are, as they should be, the same for the AC and B bond chains. They differ by 5% or less from experimental bond distances and improve upon the bond distances for pyrophyllite calculated by Cygan et al. (2004). These encouraging results indicate that the geometry optimizations meet accepted standards and, importantly, that the optimized crystal structure was not too “thin“ (Bickmore et al., 2003), i.e. the crystal structure was unperturbed by the creation of edge surfaces and large enough to isolate them from one another.
Bond distances in the two candidate edge surface structures, geometry-optimized in vacuo, are compared with ab initio calculations at high water coverage (Churakov, 2006; Table 1b), again showing agreement within a few percent. The O–H bond distances are the average for each edge face (Table 1b). The Si1–OH and Si2–OH bond lengths of the B edge are also averaged (Table 1b). These B edge Si–OH sites are equivalent because the MD equilibration was sufficiently long that the final geometry-optimized structure exhibited the dynamic disorder in the positions of >Al-OH2 and >Al-OH groups reported in ab initio simulations (Churakov, 2007). The Si–OH bond lengths are somewhat shorter, while the Al–OH2, Al2–OH, and Al–OH bond lengths are somewhat longer than those predicted using ab initio simulation. The average O–H bond length on both edges is underestimated relative to the bond lengths predicted by means of ab initio simulation. The proposed edge surface models (Figures 1e,f) assume that redistribution of bond strengths creates a Lewis acid site (Sposito, 2008) at the exposed Al. The expected bond relaxation occurs for >Al-OH2, while bond contractions occur at either the transition O on the AC edge surface or for >Al-OH on the B edge surface, relative to equivalent bonds in the bulk crystal structure. The contracted length of the transition >(SiO) –Al bond is somewhat longer than that predicted by ab initio simulation, whereas the contracted length of the transition >Si– (OAl) bond is shorter than both the ab initio prediction and the bulk Si–O(apical) bond. Relative to the bulk crystal structure, edge surface Si–OH and Al–OH2 bond lengths are relaxed (i.e. longer) and the transition O and aluminol sites which originate from the inner OH group are contracted.
Surface energies (equations 2–4) were calculated for three pore widths (6, 30, and 60 Å; Table 2). Calculation as a function of pore width is necessary to create a basis for comparison with surface excess energies (Churakov, 2006) and to confirm that the simulation cell configuration is sufficient to isolate opposing edge surfaces from one another. For the basal plane and the B edge surface, an increase in the surface energies occurred with increase in pore width from 6 to 30 Å. For the AC edge, the surface free energy, , did not change with an increase in pore width from 6 to 30 Å, then declined with the pore width increase from 30 to 60 Å. The combined surface free energy and stress, , of the AC edge surface increased monotonically with increasing pore width. For both edge surfaces, the surface excess energy, , increased slightly with an expansion of the pore width from 6 to 30 Å. The minimal changes in the surface energies of the basal plane and edge surfaces with the pore width increase from 30 to 60 Å indicate that the system with a 60 Å pore sufficiently isolates the two opposing edge surfaces.
The basal plane possessed the smallest surface free energy and stress, which is of course consistent with the layer type crystal morphology of 2:1 clay minerals. The AC edge surface energy was significantly lower than that of the B edge surface. The morphology of a well crystallized pyrophyllite particle predicted using the values of and the Wulff construction (Wulff, 1901) agreed with the observed morphology of pyrophyllite crystals: the pyrophyllite crystal structure is a parallelepiped defined by the basal planes and the AC edge surfaces (Figure 3).
The surface excess energy, , for the basal plane was not calculated because all O atoms in the basal plane are valence-saturated and the pyrophyllite interlayer does not hydrate. The values of the for the AC and B edge surfaces were always negative (Table 2), indicating that chemisorption of water molecules, hydroxide ions, and protons and the low surface energies of clean pyrophyllite were sufficient to lead to negative values which lie between the values for pyrophyllite calculated by ab initio (Churakov, 2006) and electrostatic (Bleam et al., 1993) methods.
Solvent accessible surface O
Plots of the edge site-water O RDFs for the potentially reactive O sites on the AC and B edge surfaces provide information concerning the first water shells at these sites (Figure 4). The characteristic properties, rpeak, the location of the first (or second) RDF peak for an α–β atom pair, CN, and ρ, for these RDFs are provided (Table 3).
A note to assist in the interpretation of the RDFs is necessary. By assigning the singly coordinated O bonded to an octahedral cation site as a structural water molecule with SPC/E atom types, a bias in the calculation of the number of coordinated water molecules is introduced. The bias varies depending on the distance from the site of interest to the structural water molecule and the edge surface geometry. For the AC edge, the >Si2-OH and >Si1Al-O sites possess a CN bias of ~1.0 due to their proximity to the structural water molecule. The >Al2-OH site of the AC edge includes a CN bias of ~2.0, because the structural water molecules of two adjacent Al octahedra of the edge bond chain flank this >Al2-OH site at the center of the dioctahedral vacancy. The structural water molecule of the AC edge is freely exchanged with the bulk water phase. Thus, the site-water O RDF for the >Al-OH2 site must be inferred from the second shell of the Al-O (water) RDF (Figure 4b). The >Al-OH2-O (water) CN inferred from the second Al-O (water) shell on the AC edge surface has a bias of 0.84 (Table 3). Only the O of the >Si1-O site on the AC edge surface is unaffected by the bias associated with the structural water molecule, because the distance between the >Si1-O site and the >Al-OH2 site is near 4.7 Å.
With the water O CN bias due to the >Al-OH2 site taken into account, sites >Si2-OH and >Al-OH2 on the AC edge surface have ~2.0 water O within their first shell; sites >Al2-OH and >Si1Al-O have ~1.0 first shell water O. The transition O possesses the most compact first water shell and differs from the other edge O because the transition O functions only as a hydrogen bond acceptor, whereas the other protonated edge O can be both acceptors and donors. The >Si1-OH site has the most relaxed first shell relative to the other sites, with the largest ρ (4.57 Å) and water O CN (5.16). These characteristics are consistent with the >Si1-OH site being that most exposed to the bulk aqueous phase.
The edge O-water O CNs for the B edge sites all possessed a bias due to the structural water molecule. The dynamic disorder of the >Al-OH and >Al-OH2 sites results in all B edge sites being proximal to the >Al-OH2 site. The edge O-water O CN bias for >Si1-OH and >Si2-OH sites is between 1.0 and 2.0 and varies slightly with the distribution of the structural water molecules. The water O CN bias at the B edge >Al-OH site is ~1.0. Like the AC edge, the site-water O RDF for the B edge >Al-OH2 site must be inferred from the second shell of the Al-O (water) RDF. The B edge >Al-OH2 site possesses a CN bias of 0.97.
Unlike the AC edge surface, the silanol sites on the B edge possess similar edge site-water O RDFs (Figure 4d). With the water O CN bias taken into account, the silanol sites possess 1.0 to 2.0 water O within their first shells; sites >Al-OH and >Al-OH2 have ~2.0 first shell water O. The presence of water O in excess of the expected CN bias for the silanol sites indicates that these silanol sites are probably hydrogen bond donors and acceptors with respect to bulk water. The two water O in the first shell of the >Al-OH site indicate the potential for hydrogen bond donor and acceptor interactions.
A negative has not been reported in previous ab initio simulations of the pyrophyllite edge surface or in classical mechanical atomistic simulations of the disordered pyrophyllite edge surfaces. The negative values reported, however, do not necessarily indicate that dissolution at these surfaces would be facile. The physical reality of negative surface energies and the conditions that may give rise to them, as well as their prediction via ab initio calculation, have often been discussed (Overbeek, 1978; Stol and DeBruyn, 1980; Łodziana et al., 2004; Lin et al., 2006). Negative surface energies are also a component of a simple thermodynamic model for the stabilization of nanosheets relative to the bulk mineral (Lin et al., 2006). In this latter model, layered nanoparticles with negative inter-facial energies that arise from strong adsorption can be stabilized at small thicknesses if a sufficiently positive contribution to the surface free energy exists. The already low surface free energy, strong adsorption of water molecules, and nanosheet morphology of 2:1 phyllosilicates are thus consistent with the necessary conditions that give rise to negative interfacial energies for a thermodynamically stable particle.
The term “negative surface energy” has generated some confusion in the past, but Mathur et al. (2005) provided an excellent clarification of the term by noting the distinction between surface energies of single- and multi-component systems. In the single-component system where the solid is in equilibrium with a vapor of itself, all pyrophyllite faces possessed positive values. The dissociative chemisorption of water to the pyrophyllite edge introduces an additional component to the system in the form of a surface excess of water. The exothermic process of water adsorption decreases the surface energy and produces a negative . The calculated is sensitive to both the bulk thickness of the surface slab and the manner of quantifying the energy for the chemisorption of water to the unsaturated surface. A negative only emerged for bulk thicknesses >~1.0 nm. Previous simulations of the pyrophyllite edge surface that reported positive possessed slab thicknesses of <1.0 nm (Churakov, 2006) or employed disordered edge models with bulk thicknesses that ranged from 0.8 to 2.0 nm (Martins et al., 2014).
In addition to the model thickness, the is sensitive to the value of . In molecular simulations, the relative energy differences, and not absolute potential energies, are most meaningful in evaluating the performance of a potential model. Thus, efforts were made to calculate all accessible and relevant energies directly with the CLAYFF forcefield. Other calculations of the energy of dissociative chemisorption of water to the pyrophyllite surface in potential models have been described somewhat opaquely (Martins et al., 2014) with reference to previous simulations of zeolites (Gren et al., 2010). These latter simulations used reaction and lattice energies from ab initio results to calculate the energy for the dissociative transfer of a water molecule to the lattice oxygen positions of the Si and Al edge sites. In contrast to the present study, the surface reaction at the aluminum site was modeled as a dissociative hydroxylation reaction that would produce >Al-OH and >Si1Al-OH surface groups (Martins et al., 2014). An estimated value of from the ab initio simulations of the hydroxylation reactions at the surface Al and Si is −381.07 kcal/mol (Gren et al., 2010). If the reaction of the H2O molecule at the edge Al is considered as a condensation reaction (i.e. the gas phase energy of water less the heat of vaporization of liquid water), then the revised estimate of would be −497.52 kcal/mol. The CLAYFF estimate of is bracketed by these two density functional theory (DFT)-derived estimates of and both estimates result in a negative for the AC edge; the B edge becomes positive for the lower estimate, but still remains in the range of thermal vibration at <0.90 kcal/(mol Å2). Thus, the reported negative values calculated with CLAYFF are quite reasonable.
The short MD simulations were used to assess which AC edge surface O were sufficiently near the interface to interact with the aqueous phase. The presence of an ordered water O shell at the >Al2-OH site (Figure 4b) with an average distance similar to the other surface O sites suggests that this inner OH site on the AC edge surface is indeed solvent-accessible. This result is counter to previous ab initio simulations of the edge that have neglected or dismissed this site as solvent-inaccessible (Bickmore et al., 2003; Churakov, 2007). Perhaps the longstanding neglect of this site has contributed to the unexamined assumption that this structural hydroxyl is solvent-inaccessible. Yet, previous ab initio simulations have identified transient hydrogen bonds between this exposed structural OH and the solution phase (Churakov, 2007) and predicted an acid-base reactivity intermediate to the two different silanol edge groups (Churakov, 2006). Experimental evidence for the solvent accessibility of this site comes from the hydrothermal synthesis of deuterated specimens of natural pyrophyllite. The substitution reaction of H → D yields only 5% in natural specimens (Russell et al., 1970). This low yield for pyrophyllite is similar to that for talc, the non-swelling trioctahedral analog of pyrophyllite. Deuterated talc specimens yield 2–10% substitution with exchange occurring only at the talc edge surfaces and no exchange of the H from the hexagonal ring of the talc basal surfaces (Ferrage et al., 2003). Thus the results of ab initio simulation, experiment, and the classical mechanical simulations suggest that the exposed structural hydroxyl is solvent-accessible.
The role of the O atoms in acid-base reactions and surface complexation cannot be evaluated conclusively using the short MD trajectories, but the extent of hydration and the position of a site relative to the mineral-solvent interface can be used to infer the likelihood that a site will participate in adsorption phenomena. The accessibility of the structural OH group to bulk water suggests strongly that this site should be included in models of adsorption by 2:1 clay mineral edge surfaces (Bickmore et al., 2003; Tournassat et al., 2004; Bourg et al., 2007; Dähn et al., 2011; Churakov and Dähn, 2012; Kremleva et al., 2012). The >Si-OH sites were solvent-accessible, but not equivalent, on the AC edge surface. Thus all five edge O sites of the AC edge are deemed solvent accessible. Note that the most parsimonious model of 2:1 clay mineral edge surface acid-base chemistry, constrained by titration data and the crystal structure and employing the best available pKa estimates, concluded that the edge surface should have five dissociable protons (Bourg et al., 2007).
The absence of the B edge surface from the pyrophyllite nanoparticle implies that this edge surface is not energetically favorable at the p.z.n.p.c. This absence does not, however, preclude the possible existence of the B edge in other 2:1 dioctahedral phyllosilicate minerals. The B edge surface is distinct from the AC edge surface, in that the surface O atoms occupy one plane, the >Si-OH groups are equivalent, and the edge Al possesses two singly-coordinated O. These differences result in a B edge surface that is dynamically disordered relative to the AC edge surface and one that probably possesses edge sites that are distinct from the AC edge. An intermediate state of the exchange of the structural water and the aluminum hydroxyl group is evident when the B edge surface >Al-OH group is positioned between the O atoms of >Al-OH2 and >Al-OH sites (Figure 2b). The transient exchange of structural water molecules creates an edge Al with five-fold O coordination, in agreement with previous ab initio simulations of the pyrophyllite B edge surface (Liu et al., 2012; Tazi et al., 2012).
These simulations have demonstrated that the potential-based CLAYFF forcefield can reproduce the structures and energetics of the pyrophyllite edge surface as predicted by ab initio simulations. These results pave the way for future atomistic simulations that explore surface phenomena (e.g. ion adsorption, diffusive transport in clay micro- and meso-pores, stabilization of soil organic matter, and ternary surface complexes) at greater length-and time-scales than are currently possible with more computationally expensive methods.
The geometry optimizations indicate that the AC and B edge surfaces (Figures 1e,f), modeled with CLAYFF, can provide reliable predictions of pristine edge surface bond lengths, surface energies, and nanoparticle morphology. The atomistic model of the bulk pyrophyllite crystal structure predicted bond lengths that agreed within ~5% with experimental and ab initio simulation results (Lee and Guggenheim, 1981; Refson et al., 2003; Churakov, 2006), and the predicted relaxation of the AC and B edge surface bonds agreed with results from an ab initio simulation (Churakov, 2006). The predicted surface free energy and stress, , increased in the order: basal plane < AC edge surface < B edge surface, while the predicted crystal morphology incorporating values in the Wulff construction reproduced successfully the observed morphology of pyrophyllite crystals. Negative values of the surface excess energy, , for the AC and B edge surfaces resulting from strong adsorption of water by under-coordinated edge Si and Al provide evidence in support of the physical reality of negative surface energies and a simple thermodynamic model for the stabilization of other layer-type nanoparticulate minerals (Lin et al., 2006). Each AC edge TOT unit possesses five solvent-accessible O. The predicted solvent accessibility of exposed structural OH groups (>Al2-OH), supported by proton titration data, differs from conclusions drawn in previous studies (Bickmore et al., 2003; Churakov, 2007).
Atomistic simulations were completed using the 576-processor computing cluster ‘dino’ at the MGCF, University of California at Berkeley, which receives grant support from the U.S. National Science Foundation (CHE-0840505). The authors thank Dr Kathleen Durkin, Director of the MGCF, for her assistance and many helpful discussions. Financial support for this research was provided through graduate student fellowships from the Kearney Foundation of Soil Science, Carolyn Meek Memorial Fund, and the James P. Bennett Agriculture Fund, University of California at Berkeley, and the National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Japan. Additional funding for the completion of this manuscript was provided to A.G.N. through the project: ‘Multidisciplinary investigation on radiocesium fate and transport for safety assessment for interim storage and disposal of heterogeneous waste,’ under the Initiatives for Atomic Energy Basic and Generic Strategic Research by the Ministry of Education, Culture, Sports, Science and Technology of Japan, and to the G.S. by the University of California at Berkeley under the tenure of a Chancellor’s Professorship.
- Received 12 November 2014.
- Revised 2 August 2015.