Computational Molecular Basis for Improved Silica Surface Complexation Models

Surface Complexation Modelling (editor) J. Lutzenkirchen (editor) 2006 Elsevier Elsevier Ltd. Ltd. All All rights rights reserved. reserved. © 2006

359 359

Chapter 13 Computational molecular basis for improved silica surface complexation models N. Sahai" and K.M. Rosso" department of Geology and Geophysics, 1215 West Dayton Street, University of Wisconsin, Madison, WI 53706 b

Chcmical Sciences Division, Pacific Northwest National Laboratory, P.O. Box 999, MS K8-96, Richland, WA 99352 1. INTRODUCTION Understanding fundamental aspects of surface reactivity in solution begins with the identification of types and numbers of surface sites, and characterization of their thermodynamic stability and kinetic reactivity. For silica polymorphs and silicates in aqueous solution, we are primarily concerned with surface site reactions involving water species such as protonation/deprotonation and cation/anion adsorption. Taken collectively, these reactions dictate the development of surface charge, potential and reactivity, and thus underpin the overall surface chemical behavior of the silica polymorphs. This, in turn, is linked to numerous processes of relevance to a wide range of fields. Examples include the effects of particulate quartz and asbestos dusts inhaled into the lungs, fracture initiation and propagation in glasses, and the synthesis of silicate-based materials with potential medical, technological and industrial applications including orthopaedic and dental prosthetic implants, silicaencapsulated targeted drug delivery, nanopatterned semiconductor, photonic and optical devices, molecular sieves, colloids and catalysts. Moreover, because silicate minerals dominate the chemical composition of rocks at the surface of the earth and the oceans, reactions at silica(te)-water interfaces control the efficiency of biogeochemical and environmental processes such as the global cycling of nutrients, the dissolution and precipitation of minerals, and the subsurface mobility of contaminants. Ever since the emergence of surface complexation models (SCMs) in the 1970's and their continued development into this century, significant efforts have been devoted to assuming or rationalizing the different types of surface sites present, and the corresponding surface protonation equilibrium constants (pA^s, where p = -logio). Attempts to predict intrinsic surface acidities (and, thus, surface charge) include correlations of experimentally obtained pA^as values with (a) experimental aqueous hydrolysis constants for the corresponding metal ions, (b) hydrolysis constants calculated for small metal-water clusters using ab initio quantum chemical or atomistic computational approaches, and (c) effective charge densities calculated using various schemes for the surface site being protonated/deprotonated. These correlation-based approaches have been successful for some metal oxides, but silica is notoriously atypical, where the silica surface is found to be much more acidic than the silicic acid monomer in solution [1, 2]. It is noteworthy that the unusual behavior of Si, a metalloid, versus metals such as Al, Ti, Cr, Mn, Fe and Zn is also reflected

360

N. Sahai and K.M. Rosso

in the solution hydrolysis behavior of these elements so that silicon forms silicic acid whereas the harder metals form oxides and hydroxides [3, 4]. Attempts have also been made to relate oxide hydrolysis (dissolution) rates as a function of pH to the magnitude of surface charge [59] (see also the chapter by Mielczarski and Pokrovsky in this book). An alternative tack to the various correlation methods is to use a first principles (ab initio) approach applied to local clusters of atoms or to an extended crystal system along with a sufficient number of water molecules. In the latter case, both long- and short-range interactions would be taken into account, and one would predict the surface site types formed on a fresh, cleavage surface followed by hydration/hydroxylation, and the corresponding protonation constants. Some first steps are currently being made along these lines (e.g., [1013], see also chapter by Bickmore et al. in this book). It is our primary intent in this chapter to relate the surface site types, acidity constants and hydrolysis pathways as predicted by various molecular models for sites identified by spectroscopic methods, with the corresponding parameters required for SCM's. We will provide a description of various computational modeling methodologies and strategics, and review model predictions for the acidity and dissolution of amorphous silica and quartz surfaces. Advances due to application of state-of-the-art Car-Parrinello Molecular Dynamics (CPMD) method are highlighted. Other modeling approaches such as linear free energy relationships (LFER's), continuum and bond-valence methods are also discussed. Finally, we draw links between the results of molecular modeling studies with the SCM's. We note that, in addition to surface site types, site-densities and site-acidities, all SCM's also require some pheonomenological model for the electric double layer (EDL) at the solid/solution interface such as the Gouy-Chapman, Stern-Graham and Triple-Layer models. These models relate surface charge and potential and may involve capacitances, interfacial dielectric constants, etc. We will not address recent advances in modeling the structure of the EDL at the oxide/water interface (see chapters by Attard, Fawcett and Ohshima in this book) nor attempt any links to phenomenological models of the EDL. 2. METHODOLOGY 2.1. Computational Thermochemistry The Gibbs' free energy of reaction in gas-phase (AG°gas) is defined as (1) where T is the Absolute Temperature (K) and the standard state enthalpy is given by Atf°gas = A£gas + A£ZPE + A// herm

(2)

in which AE is the internal energy at absolute zero temperature, AEZPE is the zero point energy (ZPE), Aftf"16"" is the thermal contribution to AH°gas consisting of the vibrational, translational and rotational energies, and the standard state entropy of the reaction (AS*gas) is the sum of vibrational, rotational and translational entropies ASVs = A1S°vib + AS<)rot + ASV n

(3)

The standard state should be defined clearly in each study, as various standard states are in common use (e.g., [14-16]). Most ab initio methods directly yield a value of £gas,; for each

Computational Molecular Basis for for Improved Silica Surface Complexation Models

361

species "f' in the reaction. In most ab initio codes, values of EZPEt, /f"16™, and 5*/ are typically computed automatically after performing a frequency calculation. The overall reaction energy is then calculated in the normal fashion as AEgas = S ViEg!iSj where the summation is over all species "i" in the reaction and v-, is the stoichiometric coefficient of the z-th species which is positive for reactants and negative for products. The ZPE and thermal contributions are often negligible compared to the internal energy so that the enthalpy is dominated by the E term. Similarly, entropic contributions are often small. In many computational studies, therefore, it is assumed that relative trends in properties and energies for each molecule and for the overall reaction can be predicted with sufficient accuracy by the simplifications A£gas ~ AHgas, ~ A / / ^ , and AG°gas« A//gas » A£gas. For reactions in aqueous solution, solvation free energies (AG°S) must be added to AGogas in order to obtain the Gibbs' free energy in solution (AG°aq). We will discuss these topics in more detail below with relevance to silica systems. 2.2. Modeling Silica Systems 2.2.1. Cluster Versus Periodic Models The simplest cluster models of silica surfaces that have been used at the ab initio level consist either of one silanol bond, HsSiOH [17, 18], or, more commonly of one and two tetrahedral units, e.g., H4SiO4 and H6S12O7 [19, 20]. Polymeric clusters with eight or more tetrahedra have also been used [21]. Cluster models turn out to be very successful for simulation of certain properties or processes in silica, which include charge transfer and various related electronic responses such as exciton formation [22], structural properties and force constants [20], and electron density distributions [23, 24], Reaction energies and activation barriers for silica surface chemistry are computed with mixed success using cluster models. For example, energy minimizations using one water and one monosilicic acid molecule suggested that a slight expansion of siloxane Si-0 bonds precedes hydrolysis [25], a prediction later substantiated using more accurate periodic slab models [26]. However, a prediction of the deprotonation enthalpy of a silanol surface group using a small cluster model was less successful, resulting in large errors on the scale 100 kJ mol~ [21]. The most obvious deletion in the cluster approach is the effect of the long-range electrostatic potential of the extended lattice, although this can usually be overcome with embedding in a point charge field. However, the elastic field of the lattice around a region of the surface, which can contribute to the local relaxation or reconstruction behavior, is absent in the cluster representation. This deficiency was blamed for overestimation of activation energy barriers for dissociative adsorption of H2 onto small Si clusters designed to mimic Si (001) surfaces [27]. Large clusters are used mostly to include multiple site types, and then progressively larger ones can include regional surface site geometry to address steric constraints. In doing so, however, it has been learned that increasing the cluster size within practical limits generally does not fully recover the effects of the elastic field of the extended lattice. Using clusters of up to five tetrahedral units, Pelmenschikov et al. [28] demonstrated that connectivity of O3Si-O-SiC>3 bridges to the surrounding solid imparts steric constraints on the relaxation of the activated complex for hydrolysis, thereby modifying the height of the activation barrier. The geometry of the activated complex is substantially different from the equilibrium geometry of isolated O3Si-O-SiC>3 clusters [29], hence reaction-induced regional strain appears to play an important role in determining the rate of hydrolysis. For amorphous silica, several hundred atoms or more may be required to properly represent the surface [30], whereas for quartz surfaces a compact unit cell can be used because of minimal structural relaxation [31]. A drawback of the periodic slab approach is the introduction of long-range

362

N. Sahai and K.M. Rosso

order for surfaces where none is expected, such as for amorphous silica. This artifact is handled by attempting to obtain convergence of computed properties with respect to increasing supercell sizes. While the choice between cluster and periodic models is casespecific, the many successes of the cluster approach point to the fact that the chemical behavior of amorphous and crystalline forms of silica is often governed by short-ranged interactions [20, 23, 29]. 2.2.2. Levels of Theory A variety of high-order, computationally expensive methods have been applied to small cluster models for accurate calculation of certain properties in the silica system. Many introductory texts to these modeling approaches are available {e.g., [32, 33]). Currently there appears to be no assessment of the accuracy of various ab initio theories specifically for the silica system. Highly accurate Moller-Plesset (MP) perturbation theory is commonly employed at the second order (MP2) and has been used up to the fourth order (MP4), as well as coupled-cluster (CC) theory. But at these levels, usually only "single point" energies (energies for a given geometry) can be obtained on very small clusters optimized at some lower level of theory [18]. The most common approach at the time of this writing is the application of density functional theory (DFT) because it provides sufficient accuracy along with computational efficiency. DFT may be implemented in both the molecular (local) orbital approach (typically using Gaussian orbital basis functions) as well in the crystalline (Bloch) orbital approach using a planewave basis set. As will be discussed briefly below, the latter intrinsically involves the use of periodic boundary conditions and has been useful for modeling the chemistry of silica surfaces. Concerning the planewave DFT approach, and particularly with respect to modeling surfaces of silica polymorphs, ab initio molecular dynamics (AIMD) is proving to be a valuable recent computational development. The scheme yielding efficient AIMD was developed by Car and Parrinello [34], and this type of AIMD is referred to as CPMD. They showed that the minimization of the total energy with respect to the total wavefunction and atomic coordinates can be solved simultaneously, and that molecular dynamics can also be performed at only a slight increase in computational expense. CPMD is highly efficient quantum mechanical molecular dynamics at the DFT level of theory. The basis functions are a set of continuously oscillating planewaves traveling throughout the simulation cell. Planewaves are the most general of possible choices for basis sets, with no need for parameterization, in contrast to Gaussian local orbital basis sets. The number of planewaves depends on the cutoff energy; higher cutoff energies yield more planewaves and a better description of the wavefunction. The scattering properties of the atomic cores are often described using pscudopotentials to reduce the number of planewaves needed, thereby reducing the computational expense. Planewave pseudopotentials may be designed to be local or non-local and also, possibly, norm-conserving. Local pseudopotentials are direction independent spherical potentials that are computationally the least expensive but are also the least accurate. Non-local pseudopotentials have a dependence on the angular momentum of the valence electrons and are in more common use because of their greater accuracy. Normconserving pseudopotentials ensure that the electron density ascribed to the core leads to a match between core and valence wave functions. The current standard is the non-local normconserving form by Kleinman and Bylander [35], or the non-local form given by Vanderbilt [36]. Details of this approach can be found in several reviews (e.g., [37]). In this chapter, AIMD can be taken as synonymous for CPMD, but the reader should be aware that other approaches to AIMD exist. In the past 10 years, AIMD provided many new insights into the

Computational Molecular Basis for for Improved Silica Surface Complexation Models

363

chemistry of silica systems, including the dynamics of dehydroxylation/silanization of silica surfaces [38-40], the hydration of silica surfaces [30], and the reconstruction and hydration of a-quartz surfaces [26, 31]. We will review the findings of these studies below. Incremental improvements to the accuracy of DFT come in the form of improved exchange-correlation functionals. Early functionals were based on the limiting behavior of the uniform electron gas in what is known as the local density approximation (LDA). More refined exchange-correlation functionals include terms involving the spatial gradient of the charge density in what is known as the generalized gradient approximation (GGA). GGA gives a better description of the exponential decay of electron density in the low density region. The need for gradient corrections to the LDA in the silica system arises for the treatment of silica surfaces. While it has been shown that gradient corrections are necessary for correctly describing a high-pressure phase transition in silica, it generally does not improve the accuracy of silica phase structural properties compared to the LDA [31]. Furthermore, Garofalini and co-workers established that the LDA yielded the correct relative stabilities for the five major crystalline polymorphs of silica, a, /J-quartz, a, /?-cristobalite, and stishovite, at room temperature [41]. Given that gradient corrections make improvements to the description of the tails of the wavefunctions, it is not surprising that they yield more accurate descriptions of hydrogen bonding in ice and water [42, 43], and more accurate dissociation energies of small molecules [44]. Therefore, gradient corrections are often implemented in theoretical studies targeting the surface chemistry of silica and quartz. The BLYP [45, 46], PW [47], and PBE [44] functionals in particular are often used because they provide a good description of the hydrogen bond [48, 49]. Becke's three parameter hybrid functional (B3LYP), incorporating a prescribed amount of Hartree-Fock exact exchange, is also in common use for calculations involving silica surfaces, especially for cluster calculations of deprotonation energies [50]. Of the studies we review below involving the use of DFT, either with local Gaussian orbital basis sets or in the Car-Parrinello planewave scheme, the vast majority incorporate gradient corrections to the LDA. 2.2.3. Example Computational Schemes - Acidity The deprotonation energies of functional groups on silica surfaces are an important target for theoretical calculations. A typical strategy for adaptation of Eqs. (l)-(3) for this purpose starts with defining the reaction AHaq <-> A~aq + H+aq

(4)

where A is the acid species, and pKa = AG°aq,DP/ 2.303RT

(5)

where AG°aq?Dp is the standard Gibbs free energy for the deprotonation (DP), R is the Ideal Gas constant and 2.303 is the conversion factor from natural to base 10 logarithm. Obtaining sufficient accuracy in the calculation of AGoaq>Dp is a long-standing problem in computational chemistry. Calculations for gas-phase deprotonation reactions historically have given better agreement with gas-phase experiments, particularly for molecular species. The problem, then, is thought to lie in large part with properly accounting for solvation effects. Thus, it is potentially advantageous to separate out artificially the solution-phase reaction into a gasphase reaction with separate solvation energy contributions according to the thermodynamic cycle

364

N. Sahai and K.M. Rosso

AG< gas

HAgas

gas

H gas

5 o

9 i

CD <

HAaq AG°

aq

where AG°S is the standard free energy of solvation calculated for the geometry optimized in the gas-phase. In this scheme, one takes the upper path using gas-phase molecular calculations and includes the solvation contributions of the separate species. The gas-phase deprotonation energy can be computed using an adaptation of Equations 1-3 [17] — A n gas — 7A.3

AS° g a s =

(6)

gas

+ AS ro , + AS A£ Z P E + A£ vi

(7) 5)RT

(8)

where the entropy and energy terms involving vibrational and rotational motions are defined as differences between the acid species A" and the parent species HA, Strans(H+) ' s m e translational entropy of the proton, V2(AnI0t + 5 ) ^ r comes from the rotational and translational energies of all species involved (RT12 per degree of freedom) and the volume work pAV = RT, and Anrot is the difference in the numbers of rotational degrees of freedom ( = -1 for linear species, = 0 otherwise). Various statistical thermodynamic partition functions are available to compute the needed thermochemical terms. Many of these are automatically computed in standard ab initio codes upon calculation of vibrational frequencies. As noted earlier, careful attention must be paid to the standard state, which in many ab initio codes is referenced to 298.15 K and 1 atm appropriate for an ideal gas. For eventual application to solution conditions, however, the conventional standard state usually involves concentration units in molarity. In this case, the correction term to convert the concentration units at 298.15 K is 7.9 kJ mor 1 [15, 16, 51]. To achieve the highest level of accuracy, typical modeling considerations apply. Molecular geometries, energies, force constants, etc. depend on the size of the basis set, the level of theory in terms of inclusion of electron correlation-exchange, and assessment of the basis set superposition error (e.g. [17, 52, 53]). The geometry and energy of a molecule converges toward a constant value as the complete basis set limit is approached. Double- and triplc-zeta basis sets with polarization and diffuse basis functions typically perform well in this regard [54]. Combining these kinds of large basis sets with treatment at the B3LYP or MP2 level is currently a manageable level of computational expense that often gives good

Computational Molecular Basis for for Improved Silica Surface Complexation Models

365

results [55, 56]. Inclusion of all the terms in Equations 6-8 can provide an accuracy of approximately 4 kj mor 1 in the estimate of AGagaS for molecular species [17]. A major assumption in the use of the thermodynamic cycle above is that solvation does not alter the structure and relative energies of the gas-phase species. This assumption is inappropriate in many cases, such as in the deprotonation of silicic acid [15], and is likely to be inappropriate for silica surfaces [10, 11, 57]. Should the assumption hold true, then the above scheme can be useful, with most of the resulting accuracy limited by the calculation of the free energies of solvation. In this regard, explicit treatment of several 'shells' of solvent molecules in AIMD simulations is feasible [58], but capturing long-range effects this way is still unpractical at the time of this writing. Most of the precedent in the current literature for treating the effects of solvation is based on making various levels of approximation using classical methods (see references in Fu et ah [59]). Among the most popular approaches is to surround a molecule (the solute) by a surface which is embedded in a continuum of uniform dielectric constant, the so-called solvent reaction field, such as in Tomasi's Polarizable Continuum Model (PCM) [60] or in the Conductor-like Screening Model (COSMO) [61]. The reaction field is described by computing the polarization charge density at the cavity surface for a given charge density distribution in the solute. One ambiguity with those approaches lies in defining the size and shape of the cavity, and in certain cases how to incorporate short-range effects such as hydrogen bonding [62]. Cavities are usually defined by a set of overlapping atom-centered spheres using a predefined set of atomic radii, or by an isodensity surface of the solute. The solute charge density and the solvent reaction field it produces are usually both iterated until self-consistency is reached. Applying this approach to a series of related molecular species has led to computed absolute free energies of deprotonation that are tightly correlated with experimental ones [51, 59, 62-65], but 'tuning' the continuum model still appears to be the modus operandus for achieving a slope equal to one, especially for non-neutral molecules. Another strategy is to focus simply on determining whether or not linear correlations exist between calculated gas-phase deprotonation energies and experimental pA^a's, avoiding the calculation of absolute pATa's in solution. Note that in taking energy differences in the scheme diagrammed above, it is possible to benefit from a cancellation of errors associated with the individual absolute energies. Given the apparent uncertainties in the calculation of absolute p^ a 's, establishing linear free energy relationships (LFER's) with experimental pATa's can provide the same utility in terms of predictive power. LFER's have been shown to exist between predicted and experimental acidities of tetrahedral oxyacids and hexaquo cations calculated at the ab initio level in the gas-phase [53, 57, 66, 67, 68]. LFER's such as these capitalize on the fact that, for a properly chosen series of molecular species, the ZPE, thermal and entropic contributions to the free energy are quantities that are either simple constants or scale proportionately with A£gas. Similar philosophies underpin trends in the acidity of terminal and bridging hydroxyls in different oxyacids, which have been described in terms of proton affinity, O-H bond lengths, the O-H stretching frequencies and force constants, charge on the hydroxyl proton, effective residual charge on the oxygen atom, the degree of underbonding of the oxygen atom (e.g., [17, 53, 54, 69]). At any given level of theory, the size of the model cluster, and the type and number of bonds in the cluster also affect the predicted acidity. This is especially relevant for calculating the acidity of solid phases or surfaces, where larger model clusters or periodic calculations on mineral slabs would account for some of the structural relaxation and the steric constraints offered by the surrounding or underlying bulk solid, that would not be accounted for in a small cluster [21, 53, 54]. The predicted acidity of the double four-membered ring (HgSi802o)

366

N. Sahai and K.M. Rosso

at 1400 kj mol ', for example, is much closer to the experimental value of 1390 kj mol ' than the 1507 kJmor 1 value calculated for the isolated silanol molecule (H3SiOH) [21]. Furthermore, the ring shaped cluster has relatively few terminal hydroxyl bonds and more bridging bonds compared to a branched silicate tetramer (HioSi4Oi3) such that steric constraints, relaxation effects and hydrogen bonding effects in the ring cluster are probably more realistic. 3. CHARACTERIZING AMORPHOUS SILICA AND QUARTZ SURFACES 3.1. Surface sites on amorphous silica and quartz 3.1.1. Equilibrium Surface Sites Precipitation of solid phases from solution or crystallization from melts at equilibrium yields lowest energy growth surfaces. In the bulk solid, most of the silicate tetrahedra are completely polymerized. That is, the Si atoms are coordinated to other silicon atoms through bridging oxygen (Ot,r) bonds. At the ideal surface and at defect sites, partially polymerized silicate tetrahedra may occur. In the following nomenclature of sites, Qn refers to a Si atom that is bonded to "n" other Si atoms through bridging oxygens. When water interacts with the surface, the water molecule may dissociate, with the hydroxide bonding to Si atoms to form silanol (Si-OH) surface sites and the proton bonding to surface oxygens thus forming silanols or to bridging oxygens. A Q3 Si yields a single or isolated Si-OH group, a Q2 Si yields a pair of geminal Si-OH groups consisting of two silanols on a single Si atom, and vicinal Si-OH groups refer to Q3 silanols on neighboring Si atoms. Three types of surface silanols are thus recognized in the nomenclature of silica surface sites, isolated, geminal, and vicinal (Fig. 1). H

H

H

/\/\y

V 0/\ (a) vicinal silanols

OH

OH OH

OH

(b) geminal sitanols

I Si

]

o

o

I I

(c) isolated silanols

Fig. 1. Redrawn from Dijkstra et al. [70].

3.1.2. Generation of Defect Sites High energy, cleavage surfaces are obtained by fracturing or grinding the solid. Evidence for the different types of defect surface sites comes from vibrational Infra-Red/Raman and 29Si Nuclear Magnetic Resonance (NMR) studies, and Electron Spin Resonance (ESR) studies of paramagnetic defects such as radicals. Fracturing creates two new surfaces, and can occur either by homolytic or by heterolytic cleavage of the Si-0 bond. We rely here on excellent reviews of the different surface sites formed in vacuum as a result of each type of bond cleavage [71, 72]. In homolytic cleavage, one electron each from the Si-0 bond is transferred to Si and to O (Fig. 2).

Computational Molecular Basis for for Improved Silica Surface Complexation Models

Si

\X/

-]-• 0

Si ^]y

>Si-orQ 3 (-)or silyl radical site liomolytic cleavage of Q 3 site >Si-O- or Q 3 (O-) c

I /\\

/\\

NX/

S'

/\\

>Si* or Q3 (+) or silyl cation site

heterolytic cleavage of Q 3 site

/w

Si -y \

i /*\\

/\\

Si ~VA

homolytic symmetric cleavage of Q 2 site

>Si-O' or Q 3 (O-) o siloxy aiiion site

^ 11

^/° V V = ^ > J ° i Si

Si

Si

/ \ . . / \

\7

/

Si

\

/

\

\/

Si

Si

/ \ / \ OH OH OH OH

>Si=O o Q2 (O) OH

OH OH OH

\ /

\

Si

.

/

Si

/

\

/

\7

\

.

\/

\Si

Si ~\/\

Si "V^N

cP,

o^o

x

s/ / \

liomolytic asyniiiietric

d^vage of Q2 site

s/ / \

V / \

>Si"orQ2(")orsilylene site with lone pair of elections

Si ^

—

v V / \ / \

>Si(0 2 )orQ 2 (0 2 )or intersite peroxide group

\ / v/ 0 •

Si

/ \

0 \

0 /

\

Si

/

lieterolytic asymmetric cleavage of Q 2 site

0

\

>

si"

V I /\ /\

,,

charged silylene and deprotonated geininal diols Q 2 (O) 2

Fig. 2. Formation of various defect sites on silica and quartz: (a) and (b), vicinal sites at Q3 Si; (c) and (d), geminal sites at Q2 Si; (e) and (f), geminal and vicinal defects at Q2 Si (redrawn from Murashov [72] and Brinker et al. [73]).

367

368

N. Sahai and K.M. Rosso

Heterolytic cleavage refers to the transfer of both electrons to the more electronegative atom, oxygen, resulting in charged surface sites where Si is positively-charged and O is negativelycharged. Further, each type of cleavage can occur symmetrically or asymmetrically, and create surface sites where the Si atoms are bonded to three or two underlying silicate tetrahedra, i.e., create Q3 or Q2 surface Si sites. Homolytic cleavage in vacuum can yield Q3 Si silyl (=Si*) and siloxy (=Si-O*) radicals each with a single unpaired electron, Q2 silanones (>Si=O) which are surface sites with a double bond between Si and O similar to ketones, silylenes that are Q 2 Si sites with a lone pair of electrons (>Si") and Q2 siladioxirane (>Si(C>2)) sites. Heterolytic cleavage in vacuum yields Q3 silyl cation (=Si+) and siloxy anion (=Si-O~) sites, and Q2 silylene cation (>Si2+) and fully deprotonated geminal diols (>Si(O~)2). The surface radicals are highly reactive and are quickly quenched by water to form geminal Q2(OH)2 silanol sites and Q3(OH) silanols that are either vicinal or isolated [72]. Regardless of whether the silanols were formed at equilibrium or defect sites, once formed, vacuum degassing of amorphous silicas at intermediate and high temperature results in dehydroxylation of silanols with the formation of siloxane bonds, yielding new silicate ring-type surface sites. The formation of some of these ring-type sites is shown schematically in Figure 3. The extent of dehydroxylation depends on the temperature of heating, the distance between silanols and the amount of strain associated with the rings formed [74]. Conversely, reaction of water with ring-sites will eventually result in hydrolysis of the siloxane bond, ring-opening, and the formation of equilibrium isolated, vicinal or geminal silanol sites. The above different types of surface sites on silicas and silicates, including cyclic structures and radicals, are not merely of academic interest, but have significant implications for "real-world" applications. For example, it has been proposed that three-ring sites present on synthetic silicate-based bioceramics constitute the active site for heterogeneous calcium phosphate nucleation leading to epitaxial apatite growth at the bioceramic surface [75-78]. The surface layer apatite is responsible for bonding the bioceramic to the existing bone or tooth when implanted in the human body, which is precisely why bioceramics are used for orthopedic and dental prosthetic devices. Results obtained from the semi-empirical quantumchemical Austin Method (AMI) indicate that hydrolysis of the three-ring provides the lowest energy pathway for water-enhanced fracture (crack) growth in amorphous silica gels [79]. Similarly, surface radicals on freshly fractured quartz particles, created by blasting and construction activities, have been implicated in the complex set of reactions that ultimately lead to the lung-disease, silicosis [80-82]. Further, the acidity, hydrogen bonding ability, hydrophilicity or hydrophobicity of isolated versus geminal versus vicinal silanols and terminal versus bridging oxygens is important in controlling the catalytic behavior of silicas and silicates such as zeolites. Therefore, it is important that the molecular-level acidities of these different surface sites ultimately be related to the acidities of equilibrium silanol surface sites assumed in surface complexation models and to the macroscopic, average acidities measured by experimental acid-base surface titrations of the oxides. 3.2. Surface Hydration/Hydroxylation Models In the following sections, we will first review studies on the chemical reactions of water including hydration, hydroxylation and dehydroxylation (condensation and densification) on the amorphous silica surface, followed by the physisorption of water at increasing levels of surface coverage on amorphous silica. Hydration and hydroxylation at the dominant, ideal, crystallographic faces of quartz, namely, the (0001), (1010 ) and (1011) faces

Computational Molecular Basis for Complexation Models Computational for Improved Silica Surface Complexation

369

will be considered next, and finally, the reaction of water at defect sites on amorphous silica and quartz. H

H \

/

\

~

Si

00°

I

H

II

u

\ /

\

Si

/

\

/

\

H

/

1 ]

Si

0°

0°

O°

0

I

I

I

I

Si

Si

I I

Si

I

I

I

I

(a)

\/

\T

Si

Si

I

I

/\.

/ \

Sk

^-Si

/V°~7\

Fig. 3. Formation of various ring sites on silica (a) two-membered (2M) and three-membered (D2) rings from condensation of vicinal sites; (b) five-membered ring from condensation of geminal sites; (c) four-membered (Dl) from paracrystalline cluster fusion (redrawn from Chuang and Maciel [88] and Brinker et a!. [73]). See section 3.2.3 of text for explanation of doubleheaded arrow in (a).

3.2.1. Experimental Evidence for Hydroxylation and Dehydroxylation (Condensation, Densification) of the Amorphous Silica Surface Hydroxylatcd silica surfaces possess silanol sites that are hydrophilic because they can hydrogen bond to water molecules as proton-donor and proton-acceptor. The average density of sites determined and estimated by a variety of methods is 4.5-4.9 sites nnT2 (see Zhuravlev [83] and references therein). The dehydroxylation and densification of silica surfaces has been examined by Infra-Red and Raman vibrational spectroscopies and by 29Si NMR spectroscopy [84-89]. Vacuum degassing of silica surfaces up to ~ 170°C removes physisorbed water. At temperatures from 170-450°C, isolated, single silanols and geminal silanols dehydroxylate to form low-strain pentasiloxane rings that are part of larger

370

N. Sahai and K.M. Rosso

octasiloxane rings. Between 200 and 400°C, vicinal Q3 silanols that are hydrogen-bonded due to their proximity undergo condensation reactions yielding cyclic Q4 Si trimers. These sites are called trisiloxane or three-ring or D2 sites (Fig. 3a). Cyclic tetramers containing four Q4 Si are formed by condensation of the remaining Q3 vicinal silanols at intermediate temperatures from 450-600°C. These defect sites are also called tetrasiloxane or four-ring or Di sites. The three-rings are strained and planar, whereas the four-rings are unstrained and puckered. The D2 sites are almost absent at temperatures below 200°C, increase in number per unit area on heating up to a maximum at about 600°C and decline in number by 1000cC. The number of Di sites is less sensitive to the temperature but are most abundant below about 500-600°C. The three and four-rings are characterized, respectively, by vibrational peaks at 490 cm"1 and 605 cm"1 in the Raman spectrum corresponding to the symmetric oxygen ring breathing modes [73, 74, 84-87]. The hydrolysis reaction of Di and D? sites with water is reversible up to about 200400 c C, and can lead to the reformation of surface silanols. Hydrolysis of each D2 site produces two surface silanols [88]. On heating to ~ 600cC, the surface is believed to be dominated by Q3 silanol sites and Q4 silanols in three-, four-, and higher-membered rings, with an estimated site density of 2.2-4.5 sites nm~ for the three-rings, corresponding to 28-58% of total surface sites [87]. The three-rings hydrolysis rate is estimated to be 75 times faster than that for amorphous silica, attributed to the greater strain release associated with hydrolysis of the three-rings. The proposed mechanism involves two steps: a slow, ratedetermining step for water adsorption where H2O acts as a hydrogen bond donor to the siloxane bond in the 3-ring, followed by a faster step involving dissociation of the water molecule [87], Interestingly, when single- and double-chain silicate minerals are reacted with water, they undergo leaching of the cations up to the first 1000 A accompanied by surface protonation, surface reconstruction and densification, ultimately resulting in the formation of D] sites as characterized by the appearance of the Raman vibrational band at 490 cm"1 [90]. At temperatures greater than ~ 500-600°C only a small number of isolated silanols remain at the surface, separated by large distances, that tends to limit hydrogen bonding. The remaining subset of vicinal, isolated silanols undergoes dehydroxylation and condensation to form a highly-strained, cyclic dimer containing two Q4 silicons as edge-sharing tetrahedra (Fig. 3a). This type of surface site is also called a disiloxane or two-membered (2M) ring and is characterized by Infra-Red vibrational bands at 888 and 907 cm"1 and are strong Lewis acids towards bases such as pyridine and trimethylamine, and can be eliminated by dissociative sorption of water, methanol and ammonia [74, 91, 92]. This behavior of the strained, cyclic dimer is in contrast to normal, unstrained siloxane bonds that are generally hydrophobic, their Si atoms are not acidic and the bridging O atom is not basic. The strained, cyclic dimer is estimated to be present at a density of ~ 0.1-0.4 sites nm~ [74, 93]. 3.2.2. Overview of Models for the Amorphous Silica Surface Various models for the amorphous silica surface have been given in the literature. Garofalini and co-workers produced theoretical silica surfaces by simulating glass 'fracturing' and relaxation with atomistic MD simulations on periodic systems [94-96]. Because of the lack of long-range structural order, the resulting surfaces are complex but concentrations of certain kinds of sites have been determined. Site densities determined by Garofalini [94] are: nonbridging oxygen (1.9 nm~2), three-coordinated Si (0.6 nm"2), three-coordinated oxygen (1.4 nm"2), two-membered rings (0.13 nm"2), three-membered rings (1.4 nm"2), and fourmembered rings (1.8 nm"2). It has long been suggested in the literature that specific crystal

Computational Molecular Basis for for Improved Silica Surface Complexation Models

371

faces of yS-cristobalite and/or tridymite are good representations for patches or segments of atoms (on the order of about tens to hundreds of atoms) at the surface of various amorphous silica types including silica derived by aerosol and sol-gel synthesis methods, fumed silicas, and even silica glass and vitreous silica ([88, 97, 98] and references therein). In the idealized model, the (111) face of/?-cristobalite represents segments that contain isolated, terminal silanols and the (100) face represents local areas bearing geminal silanols. Before presenting the details of the model, we note that this idealized structure may be distorted on a real amorphous silica surface due to thermodynamic or kinetic reasons so that a range of hydrogen bond strengths may be expected. For the purposes of modeling, oligomeric silicate ring clusters have been suggested as practical models for amorphous silica because they can be synthesized to possess isolated terminal silanols, geminal silanols and vicinal silanols [70, 99], and the acidity and reactivity of these sites to water, methanol, and amines can be tested individually. The oligomeric ring clusters can also be used in ab initio cluster calculations of surface structure, acidity and reactivity because of their relatively small size {e.g., [53, 100]). 3.2.3. Modeling (111) fi-Cristobalite-like Patches on Amorphous Silica Chuang and Maciel [88] used 29Si Cross-Polarization and Direct Polarization Magic Angle Spinning Nuclear Magnetic Resonance (CP-MAS and DP-MAS NMR) studies to characterize silica surface sites in terms of the idealized yS-cristobalite model. These authors showed that the isolated, Q silanols on the (lll)-like patches cannot hydrogen bond with each other because the O-O distance between adjacent silanols (see dashed, double-headed arrow in Fig. 3a) at 5.0 A is too far compared to the 3.3 A required for hydrogen bonding. Each isolated silanol may, however, hydrogen bond to water molecules. Dehydroxylation of two isolated sites near each other would yield the strained D2 sites. Furthermore, geminal Q2 silanols belonging to the same silicon also cannot form hydrogen bonds with each other, even though they are close, because of their mutually unfavorable orientation. The inner pair of silanols from two adjacent geminal sites can hydrogen bond to each other because they are only 2.3 A apart and they have the appropriate orientation. Hydrolysis of a siloxane bond between two silicon atoms would yield identical hydrogen-bonded geminal silanols. Dehydroxylation of these Q2 geminal silanols would yield a low-strain pentasiloxane ring as part of a larger octasiloxane ring [88] (Fig. 3b). Hydroxylation, dehydroxylation and reactivity towards ammonia of the (111) and (100) faces of /3-cristobalite were studied by AIMD simulations using the Car-Parrinello scheme [38-40]. The idealized geometry of the (111) face yielded 4.55 single silanols mrf2 [38], close to the experimental value of 4.9 sites nm~ for amorphous silica [88]. Consistent with the geometrical arguments of Chuang and Maciel [88], it was found that the isolated silanols on (111) surface are not hydrogen-bonded. Dehydroxylation of the Q3 isolated silanols on the (111) face resulted in formation of 2M rings [38] with a condensation energy of 127 kJ mop' which is higher than the experimental value of 81 kJ mol"1 (Fig. 3a). The 2M rings are strained so a higher energy may be expected compared to the experimental system where both low-strain five- and fourmembered rings as well as high-strain three- and two-membered rings are present. A calculated site density of 0.5 sites nnT was obtained for the two-membered rings close to the experimental range of 0.1-0.4 sites nm~ , and five-and seven-membered rings were created just below the 2M ring consistent with Chuang and Maciel's arguments [88]. Vibrational peaks of 849 and 856 cm"1 were obtained for the two-membered ring, which if corrected

372

N. Sahai and K.M. Rosso

upwards by 20 cm ' for computational error, compare favorably to the experimental spectra at 888 and 907 cm"1 [38]. The reactivity of the 2M rings on the (111) face with water and ammonia was also calculated by AIMD in the Car-Parrinello scheme [40]. Two different hydroxylation pathways were investigated (Fig. 4a, b). In the first, the water molecule adsorbed at the silicon atom in the ring where the silicon acted as an electrophile (pathway "A"). The proton adsorbed at the bridging oxygen due to its basic nature in the second, "B" pathway. The activation energy for pathway B was 30.89 kJ mol"1, 75.25 kJ mor 1 lower than for pathway A suggesting that the chemisorption and hydroxylation reaction was driven by the basicity of the oxygen. A preceding step involving the physisorption of water at the ring occurred at the Si atom, along pathway A, with an adsorption energy of 10.61 kJmor 1 . The calculated hydroxylation and activation energies of 164.0 and 30.87 kJ moF 1 obtained by the AIMD approach are close to the earlier values of 135.07 and 27.01 kJ mol"' obtained for the minimal Si(OH)4 cluster with the DFT-BLYP approach using Gaussian basis sets [101]. Thus, it appears that, at least for the hydroxylation reaction, the results are not greatly dependent on accurate long-distance models of the amorphous silica surface. The physisorption of ammonia at the 2M rings on the (111) face was favored by a deep well of 43.42 kJ mol"1 where sorption occurred with nitrogen donating its electron pair to the Si atom [40]. As with water, the chemisorption process can occur along two pathways (Fig. 4c, d), and the lower energy pathway involved sorption of the proton from the ammonia at the bridging oxygen site ultimately resulting in ring opening with the formation of Q3 Si-OH and Q3 Si-NH2 sites (Fig. 4d). n

\sf / ° X si'/ / X ^ X 0 o""

\

OH

Oil

-CKr

> / s\r°N Q/ \

b (

C (

>

/

X O

/

X

\

OH

NH 2

Fig. 4. Model for (a, b) hydroxylation and (c, d) ammonium adsorption pathways on 2M sites (modified from Masini and Bernasconi [40]).

Computational Molecular Basis for Complexation Models Computational for Improved Silica Surface Complexation

373 373

The nature of water-2M ring interaction has also been studied using QM/MM dynamics simulations [30]. A metastable, five-coordinated Si atom was found before the water molecule dissociated, with three Si-O bonds from 1.81 (Si-Orjng), 1.82 (Si-Orjng) and 1.84 A (Si-OH2) and two Si-0 bonds of 1.64 and 1.67 A (Fig. 5). Different results were found for an isolated 31-atom cluster without the bulk silica matrix, and for an isolated 2M ring (H4Si2O4). The Si-OH2 distance was 2.56 A in the 31-atom cluster, and water did not bond to the isolated 2M ring. These results were interpreted as suggesting that the long-range interactions of the two-ring and water with the embedding silica matrix are important in predicting the correct structures and energies. On the other hand, as mentioned above, very similar results were obtained for physisorption of water at two-membered ring sites on (111) yff-cristobalite using AIMD [40], and for water interacting with a silicic acid monomer cluster [101]. Therefore, the discrepancy between the periodic and cluster calculations observed in Du et al.'s study [30] could also reflect problems in the QM/MM model, such as in the treatment of the link atoms at interface between the QM and MM regions, or in the forcefields in the classical region of the simulation.

O,H

OjH

Fig. 5. Interaction of (a) a water monomer and (b) two water molecules with the 2M ring, which is part of a 31 atom cluster treated quantum mechanically, and embedded in an amorphous silica matrix treated by classical MD (representing results of Du et al. [30]).

The dissociation of the water molecule proceeded with the proton from water attaching to the bridging oxygen of the 2M ring and the calculated energy barrier was 38.59 kJ mol"1 [30]. This pathway of initial physisorption of water at the Si atom followed by proton sorption at the Ormg is similar to pathways A and B of Masini and Bernasconi [40], and the calculated energies from the two studies are in very close agreement. The addition of a second water molecule changed the reaction pathway significantly such that the energy barrier disappeared. The first water molecule transferred a proton to the second forming hydroxide (OH) and hydronium (H3O) ions. The OH" bonded to one Si atom of the 2M ring while the hydronium transferred a proton to the bridging oxygen of the ring, Ormg. Thus, in the subsequent reaction, the first water molecule dissociated, broke the silaxone bond forming two surface silanols, while the second one remained intact and within H-bonding distance of one of the two newly-formed silanols [30].

374

N. N. Sahai Sahaiand andK.M. K.M.Rosso Rosso

3.2.4. Modeling (100) /3-Cristobalite-like Patches on Amorphous Silica Dehydroxylation of the (100) faces of/?-cristobalite has been simulated using AIMD in the Car-Parrinello scheme, and compared with results for the (111) face [39-40]. Neighboring geminal silanols on the (100) face were found to be hydrogen-bonded, unlike the isolated silanols on the (111) face [38]. Geminal sites on the (100) face where one silanol accepted a hydrogen bond and the other silanol donated a hydrogen bond were more stable than geminal sites where both silanols accepted hydrogen bonds. The difference in energy between the two types of surfaces was 7.2 kJ moP 1 per geminal silanol, comparable to the a-/? transition of cristobalite [39]. Because of the small difference in energy, both types of geminal silanol configurations were expected at the (100) face of /?-cristobalite. For vicinal silanols on the (100) ^-cristobalite face, dehydroxylation resulted in the formation of five-membered rings {e.g. Fig. 3b) for fully relaxed surfaces using AIMD in the Car-Parrinello scheme [39]. The calculated dehydroxylation energy was 42 kJ mol"1 which is about half of the experimentally observed value 8 1 k J m o r ' at 350-650cC. This was rationalized by noting that the five-membered rings are low-strain structures compared to the experimental system where the observed energy reflects the average of high-strain and lowstrain surfaces such as two-, three-, four-membered rings. The calculated value was also much lower than the values obtained in previous Hartree-Fock (HF) and semi-empirical calculations that were performed with partially relaxed systems, but was similar to a recent DFT calculation that obtained 50 kJ mol~' on a fully-relaxed surface. The activation energy was calculated at 313 kJ moF 1 as the difference between the energy of the initial and the transition states. The experimentally obtained values are 80-200 kJ mol"1 at high silanol surface coverage and 200-300 kJ moF 1 at lower silanol surface coverage. The discrepancy was attributed to the use of a small supercell, which resulted in an infinite row of geminal pairs in the calculation rather than a simple geminal pair. 3.2.5. Associative (Molecular) Water Adsorption on Amorphous Silica The interaction of increasing amounts of water coverage on the (100) and (111) surfaces of ^-cristobalite was examined by AIMD simulations [102, 103]. Water coverage was defined as the number of water molecules per number of surface hydroxyls. One and two water molecules (the water "monomer" and "dimer"), half monolayer and full monolayer cases were considered for both (100) and (111) surfaces. The full monolayer consists of four water molecules per l x l unit cell of ^-cristobalite where the m * n notation refers to lattice parameter multipliers parallel to the (100) and (111) faces. Half monolayer coverage was modeled as four water molecules in a doubled "V2 x V2" unit cell on each face, that has a total surface area equivalent to double a 1 x 1 unit cell. Geometrical and energetic results obtained for the half monolayer coverage were very similar to those for the water dimer in the single surface cell, suggesting that two cases are equivalent and that the finite size of the supercell has negligible effect on the accuracy of the calculation. Thus, the water monomer in the 1 x 1 unit cell was equivalent to 0.25 or quarter monolayer coverage. Adsorption of the water monomer on ^-cristobalite (100) was most favorable between two neighboring geminal sites where the water molecule, forming three hydrogen bonds, acted as a proton acceptor from both the geminal sites and as a proton donor to one of the pair of geminal sites (Fig. 6a). The calculated adsorption energy was 60.01 kJmor 1 . In the next most stable configuration, the water monomer adsorbed to the silanols of a single geminal site forming two H-bonds with silanols, acting as proton donor and proton acceptor (Fig. 6b). Finally, adsorption of the monomer with a single hydrogen bond to only one silanol of a

Computational Molecular Basis for for Improved Silica Surface Complexation Models

375

geminal site was the least energetically favorable at 32.79 kj mol (Fig. 6c). Compared to a free water molecule, the adsorbed water in all three configurations was significantly deformed, with one of the O-H bonds being longer and the HOH angle being enlarged slightly, as also noted for H2O adsorption on Pt (111), TiO2 (110) and some other surfaces. The water monomer could bridge between two adjacent chains of surface silanols.

H •••• '

H \

'•'

H

iy (a)

"O

v

77,

H H

H... ...0-...

H""

"'•• 0

o

OH

\

o

/••"•• H oC

0

/' o

c

L (c)

OH

\

Hi

J

V

(b) H

0'

••••0

\

SJO,

o

O^

1 (d)

Fig. 6. Associative water adsorption on (100) face of b-cristobalite: (a) geminal silanols with water monomer located above bridging oxygen, 3 H-bonds; (b) geminal silanols with water monomer located above both silanols, 2 H-bonds; (c) geminal silanols with water monomer located above only one silanol, 1 H-bond; (d) water dimmer located above the bridging oxygen, 6 H-bonds (representing results of Yang et al. [102, 103]).

Adsorption of the water dimer on /?-cristobalite (100) was also favored at the site between a pair of neighboring geminal sites, with four hydrogen bonds between the water dimer and the surface silanols and an internal hydrogen bond between the two water molecules, and one OH pointing away from the surface (Fig. 6d). The dimer adsorption energy was 72.17 kJmor 1 per water molecule, more stable than the monomer by 12.16 kJ mol~' per H2O compared to two water molecules separated by a large distance. The stabilization was attributed mainly to the extra hydrogen bonds. The hydrogen bond donor H2O was situated 0.5 A closer to the surface than the hydrogen bond acceptor. Compared to the monomer, the hydrogen bond donor of the dimer was 0.15 A closer to the surface, and the hydrogen bond acceptor was further from the surface by 0.35 A than the monomer. The 0 - 0 distance of the sorbed dimer (2.53 A) was shorter than in the gas-phase dimer (2.89 A) and shorter than the water dimer adsorbed on metals (2.7 A). The internal hydrogen bond was longer (1.04 A) compared to the gas-phase dimer (0.98 A). In brief, enhanced hydrogen bonding was observed for the water dimer as also seen on metal surfaces.

376

N. Sahai and K.M. Rosso

The dimer results are equivalent to a half monolayer coverage obtained for four water molecules in the doubled "V2 x V2" surface unit cell. Another structure was found for the half monolayer coverage case, where water formed a one-dimensional wire structure between two adjacent chains of surface silanols. This structure was less stable, however, by 9.65 kJ mol"1 per H2O than the more stable water dimer structure described above. Finally, water formed a two-dimensional, mosaic-like pattern at full monolayer coverage on the /?-cristobalite (100) face. The pattern consisted of a hydrogen-bonded network of quadrangular and octagonal rings of water molecules situated almost directly above the silanol sites. This new ice-structure was called the 2D tessellation ice structure by the original authors [102, 103]. Strong hydrogen bonds were formed within the quadrangular rings and weak ones between them. Each water molecule formed four hydrogen bonds as a double proton donor and acceptor. Each surface silanol was hydrogen-bonded to one water molecule above it. All water molecules were almost coplanar with the surface. The strongly bonded water quadrangles connected two adjacent chains of surface silanols. The adsorption energy was large, 68.69 kJ mor 1 per H2O and this configuration stable up to 300K. The 2D tessellation ice structure had a lower adhesion energy to bulk ice than to the (100) /?-cristobalite surface indicating that the structure and stability of the ice-layer was strongly dependent on the substrate. Experimental evidence for a structured layer of water at the quartz (0001) surface has been found by the use of sum frequency vibrational spectroscopy [104, 105], A different ice-structure, the Ih-\ike hexagonal structure, is formed from water on closepacked, metal substrates. In summary, one, two and four adsorbed water molecules correspond to 0.25 to 0.5 to 1.0 monolayer coverage. A significant increase in stability was obtained between 0.25 and 0.5 monolayer coverage and a smaller stabilization on going to the full monolayer on the (100) y0-cristobalite surface. The (111) /?-cristobalite surface has single, isolated silanols, and hexagonal periodicity compared the rectangular periodicity found on (100). The silanols are 5.0 A away from each other compared to the 2.3 A spacing of vicinal sites on (100). Hydrogen bonding on the (111) surface is, therefore, expected to be quite different from the (100). Water monomer adsorption was most stable when situated in the hollow site between three neighboring single silanol sites with an adsorption energy of 74.29 kJ mol"1, and the formation of three hydrogen bonds to the surface. The water molecule acts as double proton-donor and single proton-acceptor. As with the (100) surface, the hydrogen bonds between water and the (100) surface are stronger with /?-cristobalite compared to the water bonds with some other surfaces, such as Pt (111), Ag(lll), and MgO (100). For monolayer coverage, each water molecule is located in the hole between three adjacent single silanol sites, similar to the water monomer, with an adsorption energy of 67.63 kJ mol"1 per H2O. No hydrogen bonds form between H2O molecules and no tessellation ice-like structure is formed. 3.2.6. Modeling Hydration and Hydrolysis of Ideal Crystallographic Surfaces on Quartz The relative stability of the several low-index surfaces of quartz and their interactions with water and Na+ ions were studied have been studied by de Leeuw and co-workers using atomistic and AIMD simulations [106, 107]. Atomistic simulations of the vacuum-terminated surfaces showed that the (0001) surface has the lowest energy, followed by (1011), (1011), whereas the (1010) surface was the least stable. The (0001) surface was the most stable because all the surface Si atoms relaxed and became four-fold coordinated (Q4) with the formation of siloxane bridges at the surface. Both three-coordinated (Q3), and Q4 silicons

Computational Molecular Basis for for Improved Silica Surface Complexation Models

3 77 377

were found on the other faces of quartz. The Q3 silicons are expected to be the most reactive towards water. Molecular adsorption (physisorption) of water on the (0001) surface resulted in only a slight lowering of energy, compared to the other surfaces, because the surface atoms were already fully coordinated. Hydrogen bonds are formed between the H of water and surface OH groups, and between the OH of water and surface Si atoms. The adsorption energy of water on (0001) was calculated as 75.8 kj moF1, consistent with physisorption of water. On (1010), the oxygen atom from water adsorbs to the three-coordinated Si atoms, and both H atoms from another water molecule coordinate to the singly-coordinated surface oxygen, and the hydration energy is favorable by 18kJmoi~'. Water adsorbs similarly on the (1011) surface and the remaining under-coordinated oxygens form siloxane bridges with a net favorable hydration energy of 230.2 kJmoF 1 . Finally, the hydration of (1011) resulted in larger water-Si distances and a lesser extent of hydrogen bonding so that the hydration energy of 95.1 kJ mor 1 was in the physisorption range similar to (0001). Dissociative adsorption of water was much more favorable for all the four surfaces, and the relative stability of the four faces obtained was different. In order of decreasing stability, the hydroxylated surface energies were (1011) < (1011) < (1010) < (0001). The hydroxylated (0001) surface was the only surface that had geminal silanols. The silanols did not interact with each other but, rather, with the silanol of the vicmal, geminal group with an O-H distance of 2.08 A. The (1010) surface had single silanol sites that were vicinal to each other, and more closely spaced at 1.95 A so that the hydrogen bonding was more extensive. The greater degree of hydrogen bonding may explain the large hydration energies, which indicate that desorption of water from the hydrogen-bonded (0001) surface would be more difficult. The hydroxylated (1011) and (1011) had isolated silanols that did not hydrogen bond with each other. Removing isolated silanols was difficult, consistent with the large hydration energy of these surfaces. The AIMD results for the (0001) face showed that, when fully-relaxed, six-membered rings possessing siloxane bridges were formed on the surface [107]. In general agreement with the atomistic studies [106], water adsorbed associatively (molecular adsorption) when the surface was characterized by siloxane bridges with low hydration energy of- 40 kJ mor 1 , and dissociative (chemisorption) was noted when under-coordinated silicon and oxygen atoms are present on the freshly cleaved surface with a much larger favorable adsorption energy of 140.7 kJ mol"1. The sorption of Na+ was modeled using atomistic, molecular statics calculations (unlike other studies reviewed here that are based on quantum-chemical electronic structure theory) adding one Na+ ion and one OH~ ion to the cleaved surface, or one Na+ and one NaO~ per surface silicon, to model increasingly alkaline solution pH's [106]. In the first scenario, Na+ bonded with varying degrees of coordination to surface oxygens at distances of 2.15-2.45 A on (0001), 2.03-2.29 A on (1010), and 1.94-1.96 A on (1011) and (1011). All four surfaces were now more stable than the unhydrated surfaces but less stable than the purely hydroxylated surfaces. For the simulations representing higher pH, the distance between Na and surface oxygens was 1.84-2.72 A, depending on the surface [106]. 3.2.7. Hydration and Hydroxylation of Radical Defect Sites on Amorphous Silica Walsh et al. [108] used both pure QM and a hybrid QM/QM ONIOM method [110, 111] to examine the hydroxylation of the amorphous silica surface at defect sites. The ONIOM (Our

378 378

N. Sahai and K.M. Rosso

own M-layered Integrated molecular Orbital and molecular Mechanics) model treats solvation in different parts of a system with different degrees of computational accuracy from a rigorous QM method for portions (layers) of the system that are of greatest interest to hybrid QM/MM for intermediate layers and MM for layers furthest away from the part of interest. Four types of defect sites on the amorphous silica surface were considered, which they denoted as QV Q^, Q°3 and Q!3, or Qqp in general. In their notation, Q refers to a silicon atom that is bonded to "p" oxygens of which "q" oxygens are "dangling" (singlycoordinated). It is important to be aware that this is not the standard Q" notation (which we have described above), but we have retained their nomenclature for consistency with their original study. As far as possible, we will attempt to relate their notation to the more standard version as used by Murashov [72]. Thus, Walsh et al.'s [108] Q°4 denotes a Si atom that is coordinated to four other silicate tetrahedral through Obr bonds, i.e., Q4 or =Si-Ot,r-Si= in standard notation; their Q' 4 site corresponds to Murashov's [72] Q3 siloxy radical or siloxy anion sites (sSi-O' or =Si-O); Q°3 denotes the Q silyl radical or silyl cation sites (=Si' or sSi + ); and QS denotes a Q2 >Si-O~ site, for which no direct equivalent exits in Murashov's [72] scheme. The closest related site is the Q deprotonated geminal diols (>Si(O~)2). Using atomistic MD simulations, Wilson and Walsh [109] had shown previously that the most likely pairs of defect sites were Q V Q V Q V Q V Q V Q V and Q V c A (Fig. 7). These first two pairs of sites were represented, respectively, by the model clusters S12OH5 and Si2O2H5. A Si2O2H4 cluster represents the next two pairs of sites. For comparison, Xiao and Lasaga's [29] H6Si2O cluster for studying hydrolysis of the Si-O-Si bond represents a Q°4-Q04 pair of sites on an ideal, defect-free surface.

H^ I H

Fig. 7. Defect site pairs on silica (modified from Walsh et al. [108]).

Two pathways involving different transitions states were found for the Q°3-Q°4 defect pair, both leading to the same products SisHOH and SiF^OH. Thus, one three-coordinated Si remained even after hydrolysis so this type of reaction was called defect-conserving. The presence of the Q°3 defect reduced the energy barrier relative to the QVQ04 pair by 3 0 k J m o r ' . The effect of a rigid bulk matrix of silica on the barriers was examined by embedding the defect pair in a 3-membered and a four-membered ring. The results were found to be similar to those for the isolated cluster. Similar pathways were obtained for the

Computational Molecular Basis for for Improved Silica Surface Complexation Models

379 379

isolated Si2O2Hs cluster representing the Q V c A defect pair and for the cluster embedded in a three-ring. The Q V c A defect pair represented by the S12O2H4 cluster was considered equivalent to the two-membered ring. No transition state could be identified so the calculation was performed, instead, directly on a two-membered ring. A transition state with a small energy barrier was obtained. The Q V Q ^ defect pair and the two-membered ring were also placed at the edge of a cubic silica cluster and partially constrained geometries were obtained. The barrier in the constrained Q V Q ' 4 system was lower than for the constrained two-membered ring. Hydrolysis of both results in vicinal silanols with four-coordinated Si atoms as reaction products, thus consuming the defect and passivating the surface site. This reaction was called a defect-annihilating reaction. The isolated Si2OaH4 cluster and the cluster embedded in a three-ring was used to represent the Q'3-Q O 4 defect pair. There was no cleavage of a siloxane bond in this reaction, only splitting of the water molecule. A transition state with a low energy barrier was obtained for both the isolated and the embedded clusters. The reaction was deemed defect-annihilating because it produces fully-coordinated Si atoms. The Q°3-Q'4 and QVQ 1 ^ pairs were expected to hydrolyze rapidly because of the low energy barriers associated with these reactions and result in passivation of the defect sites, whereas the Q V Q V and Q V c A pairs reacted more slowly with water because of the higher energy barriers and resulted in retention of defect sites. 3.2.8. Modeling Surface Cleavage and Generation of Radical Defect Sites on Quartz The structural and electronic properties of the (0001) surface of quartz cleaved in vacuum and subsequent hydration were studied by AIMD using the Car-Parrinello method [26, 31]. Two different cleavage slices perpendicular to the [0001] direction were selected. The first, called the "cleavage surface", had one, singly-coordinated, non-bridging oxygen (siloxy radical) attached to one Q2 silicon (silyl radical) per surface unit cell, i.e. >Si'O' (modified version of Fig. 8a). The second slice was a 2 x 1 reconstruction of the cleavage surface such that a 2M ring consisting of two edge-sharing tetrahedra was formed at the surface (Fig. 8b). This surface was unstable because of electron density associated with an oxygen atom located just below the dimer. The 2 x 1 reconstruction therefore relaxed to a new surface where both oxygens of the ring formed an O2 dimer attached to a single silicon (siladioxirane), leaving the other silicon as a Q2 silylene. This surface was called the "dimer surface" and was energetically less stable than the cleaved surface (Fig. 8c). By slightly displacing the oxygen underlying the two-membered ring and allowing for more surface relaxation, a different 2><1 reconstructed surface was obtained. This new surface consisted of three-membered rings and represented a densification of the surface, so this new reconstructed surface was called the "semi-dense surface". The semidense surface was more stable than the cleaved surface because it has all Q4 Si and doubly-coordinated O atoms. Heating the dimer surface to 300 K yielded another 2 x 1 reconstruction in which all the Si were Q4, one oxygen was singly coordinated and one O was triply coordinated forming a valence alternation pair so this surface was called the "VAP surface" (Fig. 8d). Like the semidense surface, the VAP surface also had three-membered rings where the triply coordinated oxygen was bonded to a silicon in the three-ring. Based on IR and Raman spectra of amorphous silica, Lucovsky [112] estimated a high site density of 1019 sites cm"3 for these VAP defect sites per unit volume of amorphous silica.

380

N. Sahai Sahai and K.M. Rosso

S ^;si /N>9,. \

surface

Si / \

/

\Jrw\rrr J r

(c)

Fig. 8. Surface sites on the cleaved (0001) quartz surface: (a) the cleavage surface >Si-O" or Q2 (O") or siloxy radical site; (b) the hypothetical 2 x 1 reconstruction with edge-sharing tetrahedral; (c) the dimmer surface consisting of siladioxirane and Q silylene sites. Ob, indicates bridging oxygen sites. Modified from Ringnanese et al. [26, 31]). Finally, the most stable surface was obtained by heating the cleaved surface to 300 K where six-membered and three-membered rings (which do not exist in bulk quartz) were formed by linking Si and O atoms of the uppermost surface layer with those of the underlying layer. This surface had only Q4 Si and doubly-coordinated oxygens, and was called the "dense surface". The third layer of atoms from the surface were slightly displaced but layers below the third, within 5 A from the surface, returned to the bulk structure. The dense surface was more stable than the cleaved surface by 11.58 kJ mol"1. 3.2.9. Modeling Hydration and Hydrolysis of Radical Defect Sites on Quartz The cleaved surface, which had a Q2 Si and a singly-coordinated O, was hydrated by adding one H2O molecule per Si [26]. The O(H) from water bonded to the undercoordinated Si forming a Q3 silanol but the H from water did not go directly to the singly-coordinated O. Instead, it hydrogen-bonded to a neighboring under-coordinated oxygen atom. The subsequent hydration process was complex involving hydrogen bond formation with neighboring surface sites, and can be found in detail in the original study [26]. Hydrogen bond chains were formed at the surface between the hydrogen of the water molecule and a neighboring non-bridging oxygen of the surface. Thus, hydrogen-bonded, vicinal silanol sites were formed. A partially hydrated surface was represented by the "semi-hydrated" surface which had two Q3 Si(OH) groups per unit cell. Addition of a water molecule to the semihydrated surface also resulted in the formation of hydrogen bond chains between vicinal silanols. The fully hydrated surface was more stable than the cleaved surface and the dense surface, respectively, by 1360 kJ mol"1 and 443.8 kJ mol"1. Thus, the original cleaved surface was found to be extremely hydrophilic. Attempts to hydrate the dense surface, which had three- and six-membered rings, resulted in the hydration energy increasing exponentially as

Computational Molecular Basis for Improved Silica Surface Complexation Complexation Models Computational

381 381

the distance of approach of the water molecule changed from about 4.5 to 1 A, confirming the idea that siloxane linkages tend to be hydrophobic [26]. 3.3. Hydrolysis of Amorphous Silica and Quartz Surfaces Hydrolysis implies the splitting of a bridging or siloxane Si-O-Si bond by a water molecule yielding two silanol SiOH sites =Si-O-Si= + H2O = =SiOH + HO-Sk

(9)

Hydrolysis breaks up the bonding network in silicates and is, therefore, considered to be the fundamental step in dissolution of silicate and aluminosilicate glasses, ceramics and minerals. Physically, hydrolysis cannot be distinguished from some of the hydration/hydroxylation reactions involving 'dissociative adsorption' of water that we have already considered above, for example, at the 2M site [30]. Furthermore, the dehydroxylation reactions considered above are the reverse of hydrolysis reactions. In the preceding sections, however, our emphasis was on the type of surface site where the reaction occurred and the surface site formed by the hydrolysis reaction. For heuristic purposes, we consider hydrolysis reactions below in terms of reaction mechanisms and energies. Hydrolysis of silica(te)s is generally slow in the circum-neutral pH range, and is catalyzed at acidic and basic pH's. Using atomistic MD simulations on large simulation cells, Garofalini [94] showed that the dynamics of hydrolysis on amorphous silica surfaces are more complex than any single mechanism might imply. The reactions observed in many such simulations include: Water adsorption onto surface Si, temporary formation of pentacoordinate Si reaction intermediates, siloxane bond rupture, formation of non-bridging oxygens, dissociation of water molecules, formation of vicinal and geminal silanols, rupture of small rings, and water-mediated siloxane bond rupture. Despite the fact that many of these processes are concerted, reductionist-style attempts to isolate and treat individual steps at the ab initio level of theory on smaller representative models have been performed. For example, the proton- and hydroxide-promoted hydrolysis reactions of the siloxane bond were studied by ab initio cluster calculations at the HF and MP2 levels using the 3-21G(d), 6-31G(d) and 6-311G(d, p) basis sets [29, 113]. Results were comparable for clusters of different sizes containing from one up to five Si tetrahedra suggesting that protonation and hydrolysis reaction kinetics are controlled by local rather than long-range forces. Results obtained for a disilicic acid cluster, H6Si2C>7; showed that protonation (as H+ or as H3O+) of the bridging oxygen site results in lengthening (weakening) of the Si-0 bridging bonds such that the activation barrier for proton-promoted hydrolysis is reduced to 100.4 kJmol~' relative to 121.3 kJ moP 1 for the hydrolysis by water alone. A similar effect was obtained for hydrolysis of a Si-O-Al bond in an HeSiOAl cluster representing an aluminosilicate mineral [113]. The base-catalyzed hydrolysis of quartz was modeled at the MP2/6-31G level by using disilicic acid, H6Si2C>7, as the representative cluster [29]. Two different mechanisms were considered. In one, H2O attacks directly at a negatively charged, deprotonated, silanol site, and in the alternative mechanism, OH" attacks a neutrally-charged (i.e., singly protonated) silanol surface site. It was found that OH" attack at the neutral surface results in the formation of a negatively-charged surface site that is equivalent to water-attack at a negatively-charged site. In other words, the two mechanisms considered ultimately result in the formation of the same negatively-charged intermediate. Furthermore, two transition states are formed in both cases. One of the two transitions states involves a pentacoordinated Si atom that has lengthened (weaker) Si-0 bonds followed by attack of the water molecule at the

3 82 382

N. Sahai Sahai and and K.M. K.M. Rosso Rosso N.

weakened bond resulting in hydrolysis. The formation of the pentacoordinated species is the rate-controlling step with an activation barrier of 79.08 kJ moF1, similar to the H+-promoted pathway and much smaller than the barrier for hydrolysis by water alone. Larger clusters of 6 to 14 Si atoms from the (001) and (111) surfaces of /J-cristobalite were examined by using DFT energy minimizations at the B3LYP and MP2 levels of theory using 6-31G(d) and 6-31 lG(d) basis sets [28, 114]. The larger clusters were considered to account for steric constraints on the Si-O-Si bonds and their nearest neighbors, and results were compared for the disilicic acid molecule, H6Si2O7. The calculated activation energies for hydrolysis by water alone were 92.05 and 138.1 kJ mol"' for the hexameric and heptameric clusters representing the (001) and (111) surfaces compared to the 71.13 kJ mol"' obtained for the dimer. The greater activation barriers indicate that the underlying bulk solid prevents the activated complex from completely relaxing. Moreover, the larger number of Si-O-Si bridges in the (111) face compared to the (100) face provide a greater resistance to hydrolysis resulting, correspondingly, in a higher activation barrier. These results were taken to support the hypotheses that the hydrolysis of the first Si-O-Si bond is the rate-limiting step and that hydrolysis would preferentially occur at a Si-O-Si bond involving less polymerized Si atoms. Further work on a hexamer of Si showed that the calculated hydrolysis activation barrier of 121.3 kJ mol"1 was much greater than the experimentally obtained barrier of- 29.29 kJ mol"1 in the neutral pH range. The discrepancy was explained as evidence of an extremely rapid reverse reaction where the two newly-formed Si-OH groups undergo dehydroxylation and reform the Si-O-Si bond. This effect would result in a very low probability for Si-O-Si bond hydrolysis explaining why the experimentally obtained pre-exponential factor in the rateequation is so small. When this effect if considered, then the most likely pathway would be the hydrolysis of the Si-O-Si bond that is the least hindered by lattice resistance, i.e. the last bond to be hydrolyzcd, contrary to their earlier conclusion [114]. The activation barrier for such a site was calculated as 83.68 kJ mol"1, considered to be in closer agreement with experiment. Ab initio cluster calculations were performed for hydrolysis of Si-O-Si bonds in disiloxane, H6S12O and disiloxanol HjSiOSiFhOH, under conditions meant to mimic acidic and neutral pH conditions. Geometry optimizations were obtained at the HF/6-31G* level with energies at the MP4(SDTQ)/6-31G* and B3LYP/6-311+G(2d,p) levels [115]. The effect of increasing the number of one to four water molecules, and the effect of a continuum solvent using the polarizable continuum model was included self-consistently. It was found that the rate-determining step in the front-side attack of disiloxane by a water monomer or dimmer is the transfer of the proton from the nucleophilic water molecule to the leaving group (silanol). The activation barrier is lowered from 142.3 kJ mol 1 for a water monomer to 102.5 kJmol"' for a water dimer. A favorable hydrogen bond network in the reactive intermediate is possible in the presence of the water tetramer, such that the mechanism changes to back-side attack of the silicon center by the nucleophile (water). The activation barrier is only 8.37 kJ moF 1 more favorable than the front-side attack by the dimer. Pentacoordinated Si atoms are formed in the reactive intermediate for both front-side and back-side attack. The acidic hydrolysis of disiloxane was found to be more favorable by back-side attack and to involve an inversion in the configuration at silicon, with the formation of a pentacoordinated Si intermediate. The activation barrier was 92.05 kJ moF1 smaller than for the front-side attack mechanism of Xiao and Lasaga [113] although it was noted that the latter mechanism may be valid in the gas-phase. The proton-promoted hydrolysis reaction of

Computational Molecular Basis for for Improved Silica Surface Complexation Models

383 383

Cypryk and Apeloig [115] was 66.94-75.31 kJ mol ' lower than the neutral-pH hydrolysis reaction by water alone. Geometries of reactive intermediates, transition states and activation barrier for hydrolysis of the Si-O-Si bond in the disiloxanol cluster by one, two and four water molecules were found to be very similar to those in the disiloxane cluster, suggesting that the presence of the geminal silanol group had little effect on the hydrolysis rate. Further, as with the disiloxane molecule, the activation barrier decreased as the number of water molecules increased. Acid-promoted hydrolysis of disiloxanol was more favorable than hydrolysis by water alone. In the presence of the water tetramer, a stable cyclic intermediate complex was formed with cooperative proton transfer from the nucleophilic water to the siloxane oxygen such that the energy barrier was decreased dramatically to the extent that the reaction proceeds effectively without any energy barrier. The effect of including the polarizable continuum model was that energy barriers for the proton-promoted hydrolysis reaction are decreased by solvation relative to the gas-phase but that the activation barriers are relatively unaffected for hydrolysis by water alone. 3.4. Surface Protonation/Deprotonation Models The preceding discussion indicates that a wide variety of surface site types may be formed by cleavage of a fresh silica or quartz surface. These various site types are converted to vicinal (Q2) silanols or to geminal or isolated (Q3) silanols by subsequent hydroxylation, hydration, dehydroxylation, and hydrolysis. A fundamental property of these sites is their acidity/basicity or ability to deprotonate/protonate, where the relevant functional groups are a terminal silanol moeity (>Si-OH) and protonated bridging hydroxyl groups in a siloxane bond (>Si-O(H)-Si<). Silanols and siloxanes thus play an important role in silica and silicate, including zeolite, reactivity. Quantum mechanical modeling work in this area has almost exclusively been performed for molecules and cluster representations of surfaces, as opposed to periodic calculations. As was outlined in the Methodological Section above, particularly important for accurate pKa predictions is the inclusion of thermochemical contributions to the free energy and the effects of solvation. In contrast to much of the preceding discussions pertaining to periodic slab calculations, where energy quantities such as surface or hydration energies are largely based on zero-Kelvin total electronic energy differences alone, incorporation of approximate thermochemical contributions to acidity calculations for molecules is the current precedent. Calculated gas-phase values of deprotonation energies, AZS^DP often correlate well with deprotonation enthalpies (Fig. 9) [53] and, as described earlier, with experimental pATa's (e.g., [66, 116]). This was investigated at the HF and G2MP2 level for the neutral-charge oxyacids, H2O, CH3OH, HCOOH, SiH3OH, Si(OH)4, H 6 Si 2 0 7 , H3PO4, H4P2O7, H2SO4, HOC1, HCIO4, Ge(OH)4, As(OH)3, and AsO(OH)3) [53]. Linear relationships between calculated deprotonation energies and experimental pATa's were found (Figs. 10a, b). It was also found that the interaction energies of the acid anion with water in aqueous solution and interaction energies with H+ in the gas-phase are closely related (Fig. 11).

384

N. Sahai and K.M. Rosso

•

G2MP2 y a 0.785X + 5

,•

BLYP y = 0.732x + 81.74C

325

350

375

400

425

450

6-31G' HF A E d g o s (kcalmol"1)

Fig. 9. Calculated, gas-phase, deprotonation internal energy correlates with corresponding enthalpy at medium and high levels of theory (from [53]).

(a)

340

360 G2MP2 4 H d

(b)

y = -55.954 + 0.172X

380 gas

400

(kcalmor 1 )

r 2 = 0.890

6-31G' HF A E d g a s (kcalmof 1 )

Fig. 10. Gas-phase deprotonation energy and enthalpy at different levels of theory, (a) HF and (b) G2MP2, correlate with experimental solution pKa (from [53]).

Computational Molecular Molecular Basis for Improved Silica Surface Complexation Complexation Models Models Computational for Improved

385

hi other words, the solvation enthalpy for the acid anions (AHSOIVJA-) was closely proportional to AEgas.DP, which results in AEgaSiDp correlating well with experimental pA^a's. The acidity of charged acids, Si(OH)3(H2O)+, Si(OH)3O", Si(OH)2O22~ were also calculated. By correcting for the difference in solvation energy between the charged acid and its conjugate base, it was possible to include these charged species on the same A£gas,Dp-p^a correlation as the neutral acids.

340

360

380

6-31G' HF AEa

400

420

440

1

g M

(kc«l mol" )

Fig. 11. Calculated, gas-phase deprotonation energies and enthalpies correlate with experimental solution pK a 's because the gas-phase deprotonation energy correlates with differential hydration enthalpy of acid and its conjugate base anion because, i.e. water interactions with the anion and the proton scale similarly (from [53]).

y -m 3.014 * 0 324x

i 2 = 0,988

3 10

1

calc. solution iH
H+)

Fig. 12. Calculated deprotonation enthalpies corrected for hydration energies correlate with solution pK a 's (from [53]).

3 86 386

N. Sahai Sahai and andK.M. K.M. Rosso Rosso N.

The correlation predicted a value of pATa = - 5 for Si(OH)3(H2O)+. This species has not been observed experimentally, but the importance of its predicted stability and calculated pATa is important for the surface acidity of silica and quartz. Relative trends in pKR can be correlated with deprotonation enthalpy at medium levels of theory and using continuum models of solvation [53] (Fig. 12). Relative trends in acidity as described by residual charge on the oxygen, oxygen underbonding, charge of the hydroxyl proton (Figs. 13 a-c), proton affinity (PA, the negative of the enthalpy of deprotonation -A/f°gaSjDp), O-H bond lengths, and vibrational frequencies and force constants, are correctly predicted for medium to large-size basis sets at many levels of theory [53, 117]. Absolute values of PA, however, are accurate only at the higher levels of theory with larger basis sets [54]. The importance of including the various energetic, cntropic and solvation contributions to obtain accurate absolute values of AG°aq,Dp is demonstrated in modeling predictions of the stability of the pentacoordinated silicate anion, H3SiO^. The existence of pentacoordinated silicate anions as intermediate species is known in glasses, melts, some minerals, bidentate organosilicate compounds synthesized at extremely alkaline pH's (>12) and organometallic compounds [118-125] and has been proposed to be involved in the dissolution of quartz [29, 126, 127], The existence of pentacoordinated silicate anion, H5Si05~, may even exist in the gas-phase, but no experimental evidence exists for their stability in aqueous solution. The gas-phase energies of H3SiO4~-H2O versus H5SiC>5~ calculated at the HF/6-31G level were found to be too similar to reach a conclusive result of their relative stabilities [53]. Calculations on H4S1O4, its deprotonated species, and HsSiOs" were also performed at the B3LYP/6-31G** level, and at the local MP2 and HF levels with 6-31G and cc-VTZ(-f) basis sets, with and without diffuse functions [15]. Solvation contributions were calculated using the continuum, modified Born solvation approach of Ben Naim and Marcus [128] corrected to maintain consistency in standard states between the two studies denoted here as SG and BM. With these corrections, the standard solvation free energy for all solutes in the SG (2001) study [15] was AG°S,SG = AG°SJBM + 7.95 kj moP 1 . For

water, AG°SJSG was equal to 26.44 klmol" 1 . It was found that only the singly-charged tetrahedrally-coordinated anion, H3SiO4~, was stable in the gas-phase and anions with larger charges were unstable. Highly-accurate, calculated AG°aq)Dp values, within ~ 16-21 kJmoF 1 of experimental values, in the worst case, were obtained for solvated calculations of H4SiO4. The inclusion of the ZPE and thermal contributions to AG°aq,DP were thus found to be important in obtaining these accurate values. Accurate calculation of absolute acidity or hydrolysis constant values in aqueous solution must include the effects of solvation. In the SG (2001) study [15], the Gibbs' free energy of formation of the pentacoordinated anion, HsSiOs", from the deprotonated monomer, HSiCV, plus H2O was predicted to be unfavorable (AG°gas= 16.32 kJ moF1) in the gas-phase when only the internal AE and ZPE contributions were considered, but a favorable value (AG°gaS=-19.66 kJmoF 1 ) value was obtained when thermal contributions were also taken into account. But the addition of solvation free energy resulted in a final, unfavorable value of AG°gas= 22.59 kJ moF 1 for HsSiOs^ [15]. Similarly, an atomistic MD study of proton binding to periodic slab model silica surfaces indicated the importance of hydrogen bonding between surface functional groups and water molecules at the interface with solution for predicting acidity constants [57]. In that study, gas-phase pKa values for deprotonation of surface >Si-OH groups showed good agreement with experiment, whereas for protonation of >Si-OH (to form >Si-OH2+) the agreement was poor. The difference was ascribed to the lack of hydrogen bonding interaction with overlying solvent water.

Computational Molecular Basis for for Improved Silica Surface Complexation Models

387

r2 = 0.699

y = 262.346 + 119.523X

• «2° # CB 3 OH

E 400

,

/

*

(a)

HOO •

/ - . HCIO4

S^

Si{OH)4 As(OH)3 HCQOH

,3

H2SO4

•vr 1

1 2 - Z/CN

y • 278.633 + 134.311X

= 0.805

•

(b)

H) HCOOH

u. 350

Si

CH3OH

,S(OH) 4

/

H

AJO(OH>3

/

HCO t

0.6

•

monomers

o

oligomers

0.8 2-s,

(c)

0

200

400

electrostatic energy (kcalmol"1)

Fig. 13. Relative acidity trends with (a) residual charge at oxygen based on Pauling bond strengths (b) degree of underbonding at oxygen based on Brown bond valences, and (c) relation between potential at H + and electrostatic energies calculated in different ways (from [53]).

388

N. Sahai and and K.M. Rosso

Similarly, a drastic reduction in deprotonation enthalpies (-PA) of divalent metal ions is found on going from M(OH)~ to M(OH2)5(OH)~ where water molecules are added to the first coordination sphere of the metals [129]. The interaction of a water molecule with silanol was modeled at the MP2/6-31G*7/HF/6-31G* and HF/6-31G*7/HF/6-31G* levels. Water was found to interact simultaneously as a hydrogen bond donor and acceptor, where the silanol oxygen atom accepts a hydrogen bond and the silanol hydrogen atom donates a hydrogen bond [130]. The effects of hydrogen bonding, of a single water molecule around silanol and oligosiloxanediol clusters, on the acidity of the silanol group was examined at the MP4DTQ/6-31G7/HF/6-31G* level [18]. Hydrogen bond donation by the silanol group and intramolecular hydrogen bonding in the oligosiloxanediols were considered. It was found that hydrogen bond donation by the silanol group significantly increases the basicity of the silanol oxygen relative to a free silanol group, by 54-113 kJ mol"1, depending on the specific silanol or oligosiloxanediol cluster. Hydrogen bond donation results in greater electron density at the silanol oxygen thereby increasing its basic character. The effects of different models for solvation on the deprotonation energies of A13+(H2O)4, Fe3+(H2O)6 and Si(OH)4 were examined at the HF and B3LYP levels of theory [116]. Solvation energies were calculated on gas-phase optimized structures using continuum dielectric models self-consistently. Calculated deprotonation energies, taken as A£aq,DP, for each species correlated linearly with experimentally determined values of pK^. Hydrogen bond lengths were found to be shorter, from 1.6-1.7 A, when silanol was the hydrogen bond donor compared to hydrogen bond lengths of 1.8-2.1 A when water is the hydrogen bond donor. This indicates stronger hydrogen bonding in the former case, similar to the results of Pelmenschikov et al. [130]. The success of the LFER attests to the fact that for the trivalent cations, using a combination of explicit 'first shell' solvent with a continuum model for the 'second shell' and beyond works well. Extrapolation of the results for Si(OH)4 to the positively-charged, protonated species represented as [Si(OH)3H2O)]+ resulted in an absolute predicted pKa of - 2 . The actual existence of the positively charged surface species on silica is a subject of some controversy. The lack of positive surface charge at low pH's in potentiometric surface titrations is generally taken as evidence against the existence of the [Si(OH)3H2O)]+ species. There are a few recent reports, however, of the identification of the positively-charged species at extremely low pH's using X-Ray Photoelectron Spectroscopy [131, 132]. The effect of increasing the number of explicitly solvating water molecules up to 19, 10, and 9, respectively, was also considered, yielding pentacoordinated optimized geometries for A1(OH)3»9H2O, Fe(OH)3-10H2O and Si(OH)3-9(H2O)+ [116]. In general, from these collective observations we learn that the absolute predicted stability of species depends on a delicate balance of internal energy, ZPE, thermal contributions and solvation energies. Careful definition of the standard state is tacitly deduced as well.

Computational Molecular Basis for for Improved Silica Surface Complexation Models

389

4. LINKING MOLECULAR-LEVEL MODELS TO SCM's 4.1. Surface-Sites in SCM's Surface complexation models that assume all surface sites to be energetically equivalent may be classified as single-site type models. In this case, the silica and quartz surface site is always represented as a silanol group (>Si-OH). The surface silanol is amphoteric, and surface protonation reactions are usually represented as >SiOH + H+ = >SiOH2+

Kj

(10)

K^

(11)

and >SiCT + H+ = >SiOH

Sorption of metal complexes may occur at single or two surface sites resulting in monodentate or bidentate sites according to the equilibria >SiCT + M2+ = >SiOM+

(12)

and 2 >SiCT + M2+ = (>SiO)2M+

(13)

where ">" represents bonding to the underlying surface without explicitly specifying whether the Si is Q3 or Q2 polymerized. The sites in reactions 10-13 may correspond to either isolated or vicinal or geminal sites, whereas the formation of bidentate complexes may be related implicitly to vicinal sites. Protonation of a bridging oxygen in a siloxane bond (=Si-O(H)-Sis) where Si is Q4 polymerized, as in the hydrolysis of a silicate bond, cannot be represented by the silanol surface sites in the single-site type SCM's. The most widely used multi-site model is the MUlti-SIte surface Complexation (MUSIC) model [133, 134], see also chapter by Heimstra, Van Riemsdijk, et al. in this book. The energy of a site, and hence it's proton affinity, was assumed to depend on the residual charge of the oxygen atom at the site and could be related to the Pauling bond-strength of the underlying silicon atom. Two types of surface sites were considered, >Si-O~ and >Si2-O°, for which the residual charges are -1 and 0 because each silicon atom donates a bond strength of +1 to the surface O2~ anion. The latter sites would, in theory, be equivalent to a siloxane bond (=Si-O-Si=), but this type of site was considered to be inert because the estimated protonation constant turned out to be extremely low, not surprisingly, given that the site has a net zero charge. In the modified MUSIC model, the proton affinity was related to the degree of underbonding or undersaturation of valence of the oxygen atom at the surface site, where Brown bond strengths were used to calculate the bond valence of the oxygen. In this version of the MUSIC model, the (100) face of yfl-cristobalite was used as a model for amorphous silica and it was proposed that Q3 Si-OH sites in isolated pairs are present at the surface. No other types of surface sites were considered.

390

Rosso N. Sahai and K.M. Rosso

4.2. Surface Acidities in SCM's As previously discussed, the assumption that AG°DP scales as A//DP or as A£DP appears to be justified, in general, based on the fairly strong linear correlations obtained for measured solution or surface acidity constants and calculated AE or AH values for deprotonation or hydrolysis reactions [53, 57, 67, 116, 129], The pK^'s of silicic acid and silica surfaces, however, notoriously do not follow the pKa's defined by many other acids (Fig. 14). The silica surface is significantly more acidic than silicic acid in aqueous solution [2, 135] consistent with the long-noted fact that silicate oligomers are more acidic than silicic acid (Her [136], p. 214). The pA^aaq of silicic acid is ~ 9.5 compared to the surface acidity constants of silica and quartz where pKa]S is in the highly acidic regime, pA^,2S ~ 7.5, and the pH of the point of zero charge (pHPZc < 4). The values of pK^s, pKa2s and pHPZc have been estimated, respectively, as -0.7, 7.7 and 3.5 [2], -4.0, 7.9 and 2 for amorphous silica [137], and -1.2, 7.2 and 3 [2] for quartz. The negative values for the deprotonation of the positively-charged surface species, >Si-OH2+, are especially important since there are only very recent reports for the experimental identification of this site [131, 132]. Values o f - 5 and -2 have been calculated by ab initio methods for the corresponding gas-phase cluster, Si(OH)3(H2O)+ or H5SiC>4+ [53, 116]. The various versions of the MUSIC model do not distinguish between amorphous silica and quartz. The unusually low acidity of silica surfaces compared to silicic acid has been rationalized in various ways including differences in solvation behavior [2, 137, 135], bondionicity and electronegativity for silicate versus other metal-(hydr)oxides [10, 135]. A continuum, Born solvation approach which relies on the dielectric constants of the oxides, interfacial region, and the solution was used to show that silica and quartz have very low dielectric constants ~ 3-4 compared to the other oxides that have values at least an order of magnitude larger [1, 2]. Because the Born solvation energy contribution depends on the inverse of the dielectric constant, the contribution to total free energy of adsorption is very small for oxides with large values of dielectric constants. For amorphous silica and quartz, however, the inverse of the dielectric constant, and hence, solvation energy contribution is large. Further, Si(IV) is chemically 'softer' or less ionic than other metals such as Fe(III), Al(III), Ti(IV). This difference in the nature of the metals is reflected in interactions with water in both the solution phase and solid phase. Thus, Si(IV) reacts with water to form Si(OH)4 whereas the other metals form oxides or hydroxides [3]. Similarly, the interactions of water, polar solvent, with the 'softer', less ionic SiCh scales differently than with the 'harder', more ionic other oxides. The acidity of SiCh, therefore, does not lie on the same pK& versus p^ a a q correlation as the other oxides unless solvation energy contributions are accounted for [4]. Using a different approach, Bickmore et al. [10, 11] and chapter in this book) have shown that explicitly accounting for the ionicity of the metal-oxide bond, along with the undersaturation of the oxygen at the silanol site, the Si-0 bond-relaxation due to deprotonation, and a shape-factor to account for differences between octahedral versus tetrahcdral versus trigonal acids, results in an improved correlation of estimated pATagas with measured p^ a a q for a number of aqueous acids. By comparison, the modified MUSIC model treats solvation by explicit counting of hydrogen bond donation and acceptance between the oxygen at the silanol surface site and water molecules [137]. The number of hydrogen bonds, taken as 3, contributes in calculating the degree of undersaturation of the oxygen valence and, thus, the proton affinity. The increased acidity of silica surface and oligomers relative to the silicic acid monomer was examined by ab initio calculations on H4SiC>4, H6S12O7, HeSijC^, H8Si4Oi2, HioSi70i9 and HgSigC^o clusters at the HF/6-31G level in the gas-phase [53]. All oligomers

391 391

Computational Computational Molecular Molecular Basis Basis for Improved Improved Silica Silica Surface Surface Complexation Complexation Models Models

except the dimer were ring-clusters (silsesquioxanes). The correct trend of decreasing A£gasDP with increasing polymerization was predicted (Fig. 13). Different effects including the electrostatic potential at the H+ ion, electronic relaxation and geometric relaxation energies were found to contribute to the effect (Fig. 13c). In detail, the Si-O~ bond length decreased from 1.54 to 1.53 to 1.51 A in going from the Q° (H4SiO4) to the Q2 and Q3 species. 10

I

I

1

• •

1

• •

Zr(IV)

8 -

-

/

i

1

/

Km)

/

Si(IV) .

Fe(lll)/

O

Q.

6 --

-

Zr(IVj'

O CO

0)

u

-E

-

4 -

3

1

i

1.1 line +/- 1 log unit i

10

acidity in solution Fig. 14. Silica does not follow the correlation defined by other metals (modified from Sverjensky and Sahai [2]).

This resulted in larger values of the Brown bond strength (s\) at the Si atom and smaller magnitudes of oxygen underbonding 2-Sj (Fig. 13b). As the Si-O~ bond length decreases, the Si-Obr-Si bond length increases concomitantly. The dimer (Q1 Si) and the heptamer were predicted to be unusually basic compared to the other oligomers. The increased acidity was attributed to intramolecular hydrogen bonding which occurs in these two molecules because of their greater geometrical flexibility. When a limited degree of solvation was included by introducing a water molecule to the dimer, the effect of intramolecular hydrogen bonding was counteracted to some extent, such that the solvated dimer was more acidic than the gas-phase molecule [53].

392

N. Sahai and K.M. Rosso

The adsorption of background electrolyte ions at silica and quartz surfaces was modeled using analytical expressions for Born solvation and electrostatic (SE) theory [138, 139] and computationally by atomistic MD simulations [106]. At low Na+ concentrations, the simulations yielded Na-0 distances on quartz ~ 2.15-2.45 A on (0001), 2.03-2.29 A on (1010), 1.94-1.96 A on (1011) and (1011). The Na-0 distances vary from 1.96-2.45 A at very high Na concentrations representing highly basic solution pH's. Though the approaches used are very different, it is interesting to note that the effective solvated radius of the Na+ ion in solution (rej) and adsorbed at the surface (Re,j) were estimated, respectively as 1.91 and 2.42 A by Sahai and Sverjensky [139]. On the other hand, the SE model uses a constant factor of ~ 1.4 A for all cations to represent hydration at the surface, which probably does not reflect physical reality. Therefore, the close correspondence between the phenomenological SE model and the atomistic molecular statics model results is probably fortuitous but suggests, nonetheless, that the SE model is capturing at least the dominant forces involved in background electrolyte adsorption. 5. SUMMARY The acidity and reactivity of surface sites on amorphous and crystalline polymorphs of silica and other oxides control their thermodynamic stability and kinetic reactivity towards reactants in surface-controlled processes of environmental, industrial, biomedical and technological relevance. Recent advances in computational methodologies such as CPMD and increasing computer power combined with spectroscopic measurements are now making it possible to link, with an impressive degree of accuracy, the molecular-level description of these processes to phenomenological, surface complcxation models The future challenge now lies in linking mesoscale properties at the nanometer scale to phenomenological models that will afford a more intuitive understanding of the systems under consideration.

ACKNOWLEDGEMENTS NS acknowledges funding for this work from NSF Division of Earth Sciences, Geochemistry and Petrology Program, and Environmental Geochemistry and Geobiology Program through grants EAR 0208036 and EAR 0346689, and the American Chemical Society Petroleum Research Foundation for grant ACS-PRF 41777-AC2. KMR acknowledges a grant from the U.S. Department of Energy (DOE), Office of Basic Energy Sciences, Engineering and Geosciences Division, and also from the Stanford Environmental Molecular Sciences Institute (EMSI) jointly funded by the NSF and by the DOE Office of Biological and Environmental Research (OBER). Pacific Northwest National Laboratory is operated by Battelle for the DOE under Contract DE-AC06-76RLO 1830.

REFERENCES [1] [2] [3] [4]

D.A. Sverjensky, Geochim. Cosmochim. Acta, 58 (1994) 3123. D.A. Sverjensky and N. Sahai, Geochim. Cosmochim. Acta, 60 (1996) 3773. M. Henry, J.P. Jolivet and J. Livage, Struct. Bond., 77 (1992) 153. N. Sahai, Environ. Sci. Technol., 36 (2002) 445.

Computational Molecular Basis for for Improved Silica Surface Complexation Models

[5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

[25] [26] [27] [28] [29] [30] [31] [32] [33]

[34] [35] [36] [37] [38] [39] [40] [41]

393

G. Furrer and W. Stumm, Geochim. Cosmochim. Acta, 50 (1986) 1847. B.G. Zinder, G. Furrer and W. Stumm, Geochim. Cosmochim. Acta, 50 (1986) 1861. C. Guy and J. Schott, Chem. Geol., 78 (1989) 181. K.J. Knauss and T.J. Wolery, Geochim. Cosmochim. Acta, 53 (1989) 1493. W. Stumm, Chemistry of the Solid-Water Interface, John Wiley and Sons, New York, 1992, p. 173. B.R. Bickmore, C.J. Tadanier, K.M. Rosso, W.D. Monn and D.L. Eggett, Geochim. Cosmochim. Acta, 68 (2004) 2025. B.R. Bickmore, K.M. Rosso, C.J. Tadanier, E.J. Bylaska and D. Doud, Geochim. Cosmochim. Acta (in press). M. Predota, A.V. Bandura, P.T. Cummings, J.D. Kubicki, D.J. Wesolowski, A.A. Chialvo and M.L. Machesky, J. Phys. Chem. B, 108 (2004) 12049. M. Predota, Z. Zhang, P. Fenter, D.J. Wesolowski and P.T. Cummings, J. Phys. Chem. B, 108 (2004) 12061. M.D. Tissandier, K.A. Cowen, W.Y. Feng, E. Gundlach, M.H. Cohen, A.D. Earhart, J.V Coe and T.R. Tuttle, J. Phys. Chem. A, 102 (1998) 7787. J. Sefcik and W.A. Goddard, Geochim. Cosmochim. Acta, 65 (2001) 4435. R.L. Martin, P.J. Hay and L.R. Pratt, J. Phys. Chem. A, 102 (1998) 3565. J. Sauer and R. Ahlrichs, J. Chem. Phys., 93 (1990) 2575. M. Cypryk, J. Organometal. Chem., 545-546 (1997) 483. A.C. Lasaga and G.V. Gibbs, Am. J. Sci., 290 (1990) 263. G.V. Gibbs, Am. Mineral., 67 (1982) 421. J. Sauer and J.-R. Hill, Chem. Phys. Lett., 218 (1994) 333. J.Y. Song, H. Jonsson, and L.R. Corrales, Nucl. Inst. Meth. Phys. Res. B, 166 (2000) 451. K.M. Rosso, G.V. Gibbs and M.B. Boisen, Phys. Chem. Min., 26 (1999) 264. G.V. Gibbs, M.B. Boisen, L.L. Beverly and K.M. Rosso, (eds. R.T. Cygan and J.D. Kubicki), In Molecular Modeling Theory: Applications In The Geosciences, Reviews in Mineralogy and Geochemistry Vol. 42, Mineralogical Society of American, Washington D. C , 2001, pp. 345381. M. Heggie and R. Jones, Phil. Mag. Lett, 55 (1987) 47. G.-M. Ringnanese, J.-C. Charlier and X. Gonze, Phys. Chem. Chem. Phys., 6 (2004) 1920. E. Penev, P. Kratzer and M. Scheffler, J. Chem. Phys., 110 (1999) 3986. A. Pelmeschikov, J. Lesczynski and L.G.M. Pettersson, J. Phys. Chem. A, 105 (2001) 9528. Y. Xiao and A.C. Lasaga, Geochim. Cosmochim. Acta, 60 (1996) 2283. M.-H. Du, A. Kolchin and H.-P. Cheng, J. Chem. Phys., 119 (2003) 6418. G.-M. Ringnanese, A. de Vita, J.-C. Charlier, X. Gonze and R. Car, Phys. Rev. B., 61 (2000) 13250. A.R. Leach, Molecular Modelling: Principles and Applications, 2nd-edition, Prentice-Hall, Edinburgh, 1996, 744 pp. R.T. Cygan and J.D. Kubicki (eds.) Molecular Modeling Theory: Applications in the Geosciences, Reviews in Mineralogy and Geochemistry Series, volume 42, Mineralogical Society of America, Washington D.C., 2001, 531 p. R. Car and M. Parrinello, Phys. Rev. Lett, 55 (1985) 2471. L. Kleinman and D.M. Bylander, Phys. Rev. Lett, 48 (1982) 1425. D. Vanderbilt, Phys. Rev. B, 41 (1990) 7892. M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias and J.D. Joannopoulos, Rev. Mod. Phys. 64 (1992) 1045. D. Ceresoli, M. Bernasconi, S. Iarlori, M. Parrinello and E. Tosatti E. Phys. Rev. Lett., 84 (2000) 3887. S. Iarlori, D. Ceresoli, M. Bernasconi, D. Donadio and M. Parrinello, J. Phys. Chem. B, 105 (2001) 8007. P. Masini and M. Bernasconi, J. Phys. Cond. Mat., 14 (2002) 4133. F. Liu, S.H. Garofalini, D. Kingsmith and D. Vanderbilt, Phys. Rev. B, 49 (1994) 12528.

394 [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] [76] [77] [78] [79] [80] [81] [82] [83] [84] [85] [86]

N. Sahai and K.M. Rosso

C.Y. Lee, D. Vanderbilt, K. Laasonen, R. Car and M. Parrinello, Phys. Rev. B, 47 (1993) 4863. K. Laasonen, M. Sprik, M. Parrinello, and R. Car, J. Chem. Phys., 99 (1993) 9080. J.P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 77 (1996) 3865. A.D. Becke, J. Chem. Phys., 109 (1998) 2092. C.T. Lee, W.T. Yang and R.G. Parr, Phys. Rev. B, 37 (1988) 785. J.P. Perdew and Y. Wang, Phys. Rev. B, 45 (1992) 13244. M. Sprik, J. Hutter and M. Parrinello, J. Chem. Phys., 105 (1996) 1142. U. Rothlisberger and M. Parrinello, J. Chem. Phys., 106 (1997) 4658. A.D. Becke, J. Chem. Phys., 98 (1993) 1372. G.A.A. Saracino, R. Improta and V. Barone, Chem. Phys. Lett., 373 (2003) 411. A.C. Lasaga, Rev. Geophys. 30 (1992) 269. J.A. Tossell and N. Sahai, Geochim. Cosmochim. Acta, 64 (2000) 4097. J.B. Nicholas, R.E. Winans, R.J. Harrison, L.E. Iton, L.A. Curtiss and A.J. Hopfmger, J. Phys. Chem., 96 (1992) 10247. M.D. Liptak and G. C. Shields, J. Am. Chem. Soc, 123 (2001) 7314. A.M. Magill and B.F. Yates, Aust. J. Chem., 57 (2004) 1205. J.R. Rustad, E. Wasserman, A.R. Felmy and C. Wilke, J. Colloid. Interface Sci., 198 (1998) 119. T.W. Swaddle, J. Rosenqvist, P. Yu, E. Bylaska, B.L. Philiips and W.H. Casey, Science, 308 (2005) 1450. Y. Fu, L. Liu, R.C. Li, R. Liu and Q.X. Guo, J. Am. Chem. Soc, 126 (2004) 814. E. Cances, B. Mennucci and J. Tomasi, J. Chem. Phys., 107 (1997) 3032. A. Klamt, J. Phys. Chem., 99 (1995) 2224. D.M. Chipman, J. Phys. Chem. A, 106 (2002) 7413. A. Klamt, F. Eckert, M. Diedenhofen and M.E. Beck, J. Phys. Chem. A, 107(2003) 9380. G.I. Almerindo., D.W. Tondo and J.R. Pliego, J. Phys. Chem. A, 108 (2004) 166. V. Barone, R. Improta and N. Rega, Theor. Chem. Ace, 111 (2004) 237. J.R. Rustad, D.A. Dixon, K.M. Rosso and A.R. Felmy, J. Am. Chem. Soc, 121 (1999) 3234. J.R. Rustad, D.A. Dixon and A.R. Felmy, Geochim. Cosmochim. Acta, 64 (2000)1675. K.M. Rosso, J.R. Rustad and G.V. Gibbs, J. Phys. Chem. A, 106 (2002) 8133. A.T. Blades, J.S. Klassen and P. Kebarle, J. Am. Chem. Soc, 117 (1995) 10563. T.W. Dijkstra, R. Duchateau, R.A. van Santen, A. Meetsma and G.P.A. Yapp, J. Am. Chem. Soc, 124(2002)9856. V.A. Radstig, Kinet. Catal., 40 (1999) 693. V. Murashov, J. Mol. Struct., 650 (2003) 141. C.J. Brinker, D.R. Tallant, E.P. Roth and C.S. Ashley, J. Non-Cryst. Solids, 82 (1986) 117. B.C. Bunker, D.M. Haalad, K.J. Ward, T.A. Michalske, W.L. Smith, J.S. Binkley, C.F. Melius and C.A. Balfe, Surf. Sci., 210 (1989) 406. N. Sahai and J.A. Tossell, J. Phys. Chem. B., 104 (2000) 4322. N. Sahai, Geochim. Cosmochim. Acta, 67 (2003) 1017. N. Sahai, Am. J. Sci., 305 (2005) 661. N. Sahai and M. Anseau, Biomaterials, 26 (2005) 5763. J.K. West and L.L. Hench, Philosophical Magazine A, 77 (1998) 85. N.S. Dalai, X. Shi, and V. Vallyathan, J. Toxicol. Environ. Health, 29 (1990) 307. X. Shi, N.S. Dalai, X.N. Hu and V. Vallyathan, J. Toxicol. Environ. Health, 27 (1989) 435. V. Vallyathan, Environ. Health Persp., 102 (1994) 111. L.T. Zhuravlev, Langmuir, 3 (1987) 316. F.L. Galeener, Solid State Comm, 44 (1982) 1037. A.E. Geissberger and F.L. Galeener, Phys. Rev. B., 28 (1983) 3266. C.J. Brinker, R.J. Kirkpatrick, D.R. Tallant, B.C. Bunker and B.J. Montez, J. Non-Cryst. Solids,

99(1988)418. [87] C.J. Brinker, R.K. Brow, D.R. Tallant and R.J. Kirkpatrick, J. Non-Cryst. Solids, 120 (1990) 26. [88] I.-S. Chuang and G.E. Maciel, J. Phys. Chem., B 101 (1997) 3052.

Computational Molecular Basis for Improved Silica Surface Complexation Complexation Models Computational

395 395

[89] A. Burneau, B. Humbert, O. Barres, J.P. Gallas and J.C. Lavalley, (ed. H.C. Bergna) The Colloid Chemistry of Silica, American Chemical Society Advances in Chemistry Series 234, ACS, Washington D.C., 1994, Chapter 10, pp. 199-222. [90] W.H. Casey, H.R. Westrich, J.F. Banfield, G. Ferruzzi, G.W. Arnold, Nature, 366 (1993) 253. [91] B.A. Morrow and I.A. Cody, J. Phys. Chem., 80 (1976) 1995. [92] B.A. Morrow and I.A. Cody, J. Phys. Chem., 80 (1976) 1998. [93] A. Grabbe, T.A. Michalske and W.L. Smith, J. Phys. Chem., 99 (1995) 4648. [94] S.H. Garofalini, J. Non-Cryst. Solids, 120 (1990) 1. [95] B.P. Feuston and S.H. Garofalini, J. Appl. Phys., 68 (1990) 4830. [96] S.H. Garofalini, (eds. R.T. Cygan and J.D. Kubicki), In Molecular Modeling Theory: Applications In The Geosciences, Reviews in Mineralogy and Geochemistry Vol. 42, Mineralogical Society of American, Washington D. C , 2001, pp. 131-168. [97] H.C. Bergna, The Colloid Chemistry of Silica (1994). American Chemical Society Advances in Chemistry Series 234, ACS, Washington D.C., 695 pp. [98] A.P. Legrand, The Surface Properties of Silicas, Wiley, Chichester, 1998, 494 pp. [99] F.J. Feher, T.A. Budzichowski and K.J. Weller, J. Am. Chem. Soc, 111 (1989) 7288. [100] J.A. Tossell, J. Non-Cryst. Solids, 120 (1990) 13. [101] A.M. Ferrari, E. Garrone, G. Spoto, P. Ugliengo and A. Zeccchina, Surf. Sci., 323 (1995) 151. [102] J. Yang, S. Meng, L.F. Xu and E.G. Wang, Phys. Rev. Lett., 92 (2004) 146102. [103] J. Yang, S. Meng, L.F. Xu and E. G. Wang, Phys. Rev. B., 71 (2005) 035413. [104] Q. Du, E. Freysz and Y.R. Shen, Phys. Rev. Lett., 72 (1994) 238. [105] V. Ostroverkhov, G.A. Waychunas and Y.R. Shen, Chem. Phys. Lett., 386 (2004) 144. [106] N.H. De Leeuw, F.M. Higgins and S.C. Parker, J. Phys. Chem. B, 103 (1999) 1270. [107] Z. Du and N.H. de Leeuw, Surf. Sci., 554 (2004)193. [108] T.R. Walsh, M. Wilson and A.P. Sutton, J. Chem. Phys., 113 (2000) 9191. [109] M. Wilson and T.R. Walsh, J. Chem. Phys., 113 (2000) 9180. [110] S. Humbel, S. Sieber and K. Morokuma, J. Chem. Phys., 105 (1996) 1959. [ I l l ] M. Svensson, S. Humbel, R.D.J. Froese, T. Matsubara, S. Sieber and K. Morokuma, J. Phys. Chem., 100(1996) 19357. [112] G. Lucovsky, J. Non-Cryst. Solids, 35-36 (1980), 825. [113] Y. Xiao and A.C. Lasaga, Geochim. Cosmochim. Acta, 58 (1994) 5379. [114] A. Pelmeschikov, H. Strandh, L.G.M. Pettersson and J. Lesczynski, J. Phys. Chem. B, 104 (2000) 5779. [115] M. Cypryk and Y. Apeloig, Organometallics, 21 (2002) 2165. [116] J.D. Kubicki, J. Phys. Chem. A, 105 (2001) 8756. [117] J.B. Nicholas, A J . Hopfinger, F.R. Trouw and L.E. Iton, J. Am. Chem. Soc, 113(1991), 4792. [118] X. Xue, J.F. Stebbins, M. Kanzaki and R.G. Tonnes, Science, 245 (1989) 962. [119] B.T. Poe, P.F. McMillan, A.C. Angell and R.K. Sato, Chem. Geol., 96 (1992) 333. [120] R. Tacke, C. Burschka, I. Richter, B. Wagner and R. Willeke, J. Am. Chem. Soc, 122 (2000) 8480. [121] S.D. Kinrade, R.J. Hamilton, A.S. Schach and C.T.G. Knight, J. Chem. Soc. Dalton Trans., (2001)961. [122] K. Benner, P. Klufers and M. Vogt, Angew. Chem. Int. Ed., 42 (2003) 1058. [123] J.B. Lambert, G. Lu, S.R. Singer and V.M. Kolb, J. Am. Chem. Soc, 126, (2004) 9611. [124] R. Damrauer, L.W. Burgraff, L.P. Davis and M.S. Gordon, J. Am. Chem. Soc, 110 (1988) 6601. [125] R.R. Holmes, Chem. Rev., 96 (1996) 927. [126] F. Liebau, Inorg. Chim. Acta, 89 (1984) 1. [127] J.D. Kubicki, Y. Xiao and A.C. Lasaga, Geochim. Cosmochim. Acta, 57 (1993) 3847. [128] A. Ben Nairn and Y. Marcus, J. Chem. Phys., 81 (1984) 2016. [129] M. Trachtman, G.D. Markham, J.P. Glusker, P. George and C.W. Bock, Inorg. Chem., 40 (2001), 4230. [130] A.G. Pelmenschikov, G. Morosi and A. Damba, J. Phys. Chem., 96 (1992) 7422.

396 396

Sahai and K.M. Rosso Rosso N. Sahai

[131] Y. Duval Y., J.A. Mielczarski, O.S. Pokrovsky, E. Mielczarski and J J . Ehrhardt, J. Phys. Chem. B, 106(2002)2937. [132] A. Shchukarev, J. Rosenqvist and S. Sjoberg, J. Electron. Spectro. Related Phenom., 137 Sp. Iss. Si (2004) 171. [133] T. Hiemstra, W.H. Van Riemsdijk and G.H. Bolt, J. Colloid. Interface Sci., 133 (1989) 91. [134] T. Hiemstra, J.C.M. De Wit and W.H. Van Riemsdijk, J. Colloid. Interface Sci., 133 (1989) 105. [135] N. Sahai, Geochim. Cosmochim. Acta, 64 (2000) 3629. [136] R.K. Her, The Chemistry of Silica, John Wiley and Sons, New York, 1979, p. 214. [137] T. Hiemstra, P. Venema, and W.H. Van Riemsdijk, J. Colloid. Interface Sci., 184 (1996) 680. [138] N. Sahai and D.A. Sverjensky, Geochim. Cosmochim. Acta, 61 (1997) 2801. [139] N. Sahai and D.A. Sverjensky, Geochim. Cosmochim. Acta, 61 (1997) 2827.